Deep tunnel surrounding rock mechanical parameter inversion method based on three-dimensional rockburst pit contour and COA-GPR collaborative optimization algorithm

Abstract

The invention discloses a deep tunnel surrounding rock mechanical parameter inversion method based on a three-dimensional rockburst pit contour and a COA-GPR collaborative optimization algorithm, and the method comprises the steps: firstly building an FLAC3D numerical simulation model of a tunnel in a region where a rockburst pit is located, employing a Mogi-Coulomb hard rock plate crack constitutive model, and recording the total number and centroid coordinates of calculation units in a measured rockburst pit three-dimensional contour range; secondly, tunnel surrounding rock mechanical parameters serve as optimization variables, and the difference value between the total number of actually-measured damage area units and the total number of calculation units entering a yield state in the three-dimensional rockburst pit contour obtained through FLAC3D calculation is added to the total number of actually-measured units entering a yield state outside the three-dimensional rockburst pit contour and multiplied by a penalty factor to serve as an optimization objective function; then, a COA-GPR collaborative optimization algorithm is adopted, in combination with the tunnel FLAC3D numerical simulation model, an optimal solution when the optimization target function is globally minimized is searched, and reasonable tunnel surrounding rock mechanical parameters are obtained; the technical problem that mechanical parameters of deep tunnel surrounding rock are difficult to reasonably obtain through displacement back analysis due to the fact that excavation instantaneous displacement is difficult to monitor is solved.

Description

Technical Field

[0001] The present invention relates to the field of rock engineering technology, and in particular to an inversion method for mechanical parameters of surrounding rocks of deep-buried tunnels based on three-dimensional rockburst crater contours and a COA-GPR collaborative optimization algorithm. Background Art

[0002] Underground space has gradually become a new target for our exploration, utilization and development. However, as underground engineering construction continues to develop deeper, rock burst disasters have become more prominent.

[0003] Rockbursts, a typical underground engineering geological hazard, tend to occur near geologically weak surfaces such as cleavage or structural planes. The complex geological structure of these areas makes it difficult to directly and accurately measure the mechanical parameters of the surrounding rock. Rockburst intensity ratings are commonly used in engineering to distinguish the severity of rock mass damage caused by a rockburst. A key metric used in rockburst severity classification is the depth of the rockburst crater. While depth can, to a certain extent, reflect the scale of the rockburst, the surface area of ​​the rock mass destroyed by the rockburst and the volume of the crater, as multidimensional metrics, are more comprehensive and comprehensive indicators of the scale of the rockburst. The larger the surface area of ​​the rockburst and the size of the crater, the greater the volume of rock ejected, and the more severe the resulting losses. Therefore, a new method is needed: an optimized back-analysis method based on measured rockburst crater data. Specifically, feedback analysis is performed on the measured crater dimensions after a rockburst occurs to obtain reasonable numerical model parameters. This numerical model with reasonable parameters can then be used to simulate the subsequent tunnel excavation process, providing a scientific basis for the rational support design and safety and stability evaluation of deep tunnel surrounding rock.

[0004] Existing optimization back-analysis methods often use global optimization algorithms such as genetic algorithms, particle swarm optimization algorithms, evolutionary differential algorithms, and ant colony algorithms for global optimization. However, these methods require tens of thousands or even hundreds of thousands of calls to the tunnel numerical calculation model, which is time-consuming per calculation. This results in low back-analysis efficiency and is often unsuitable for optimization back-analysis of tunnel numerical model parameters due to the excessive total computational time. Combining efficient optimization algorithms, machine learning, and numerical calculations can significantly reduce the number of calls to the numerical calculation model and significantly improve back-analysis efficiency.

[0005] The Crayfish Optimization Algorithm (COA) is an emerging meta-heuristic optimization algorithm proposed by Professor Heming Jia and others in 2023. The algorithm is inspired by the foraging, summertime, and competitive behaviors of crayfish, and simulates these behaviors to find the optimal solution to the problem. Simulating natural behavior: The algorithm finds the optimal solution by simulating the foraging, summertime, and competitive behaviors of crayfish, which is intuitive and easy to understand. Temperature control: The exploration and development capabilities of the algorithm are controlled by changes in temperature, allowing the algorithm to flexibly switch between different stages. Good convergence effect: The algorithm shows good convergence effect and optimization performance in multiple test functions.

[0006] Gaussian process regression (GPR) is a machine learning method based on Bayesian theory. Bayesian regression effectively avoids the overfitting problem of other machine learning methods. By defining a function distribution, it assigns a prior probability to each possible function. The more likely the function, the greater its prior probability. GPR is more effective in solving the problem of function selection because it can evaluate the uncertainty of prediction results and adaptively obtain hyperparameters.

[0007] The basic idea behind the COA-GPR collaborative optimization algorithm is to leverage the COA algorithm's superior optimization capabilities to perform a global search (foraging) and record the location of the operator (crayfish) during the search. After a certain number of iterations (a period of searching for food), operators (other crayfish or their footprints) within a certain range of the current optimal operator (the current optimal foraging path and temperature for the crayfish) are selected from the search process (crayfish footprints) to form training samples. This is used to train a GPR local proxy model, which approximates the objective function (food) within the local range of the current optimal operator (the current optimal foraging path and temperature for the crayfish). The GPR local proxy model then predicts the optimal operator, which is then compared with the optimal operator from the COA global search, continuously updating the optimal operator until convergence (successful foraging) is achieved. Summary of the Invention

[0008] The purpose of this invention is to propose an inversion method for the mechanical parameters of surrounding rock in deep tunnels based on three-dimensional rockburst crater contours and a COA-GPR collaborative optimization algorithm. This method addresses the technical problem that during tunnel excavation under high geostress, geological structures such as faults, joints, and bedding not only alter the continuity of the rock but also significantly reduce its mechanical properties and stability, leading to sudden instability and rockburst under disturbance. Furthermore, the complexity of these geological structures leads to uneven stress distribution in the rock and localized failure modes, making it difficult for macroscopically measured mechanical parameters to accurately reflect the rock's true mechanical behavior. This method provides a scientific basis for the rational evaluation of surrounding rock stability and rational support design during the excavation of deep underground tunnel projects.

[0009] To achieve the above object, the present invention adopts the following technical solutions:

[0010] A method for inverting the mechanical parameters of surrounding rock of deep tunnels based on three-dimensional rockburst crater contours and COA-GPR collaborative optimization algorithm includes the following steps:

[0011] Step S1: Constructing a FLAC 3D numerical model of the tunnel: A FLAC 3D numerical simulation model of the tunnel in the area where the rockburst occurs is established. The FLAC 3D numerical simulation model adopts the Mogi-Coulomb hard rock slab crack constitutive model and records the total number and centroid coordinates of the calculation units within the three-dimensional contour of the measured rockburst;

[0012] Step S2: Establish the optimization objective function: Take the mechanical parameters of the tunnel surrounding rock, i.e. the mechanical parameters of the Mogi-Coulomb hard rock slab crack constitutive model, as the optimization variables

[0013] Step S3: Subtract the total number of units in the measured damage zone from the total number of units that have entered the yield state within the three-dimensional rockburst crater contour calculated by FLAC 3D, and add the total number of units that have entered the yield state outside the measured three-dimensional rockburst crater contour multiplied by the penalty factor as the optimization objective function;

[0014] Step S4: Searching for global optimal parameters: Using the crayfish-Gaussian process regression (COA-GPR) collaborative optimization algorithm, combined with the tunnel FLAC 3D numerical model, searching for the optimal solution when the optimization objective function is globally minimized, until the preset value is consistent with the measured value. At this time, the rock mass mechanical parameters in the tunnel surrounding rock mechanical parameter group are the rock mechanical parameter inversion values;

[0015] Step S5: performing forward calculations based on the FLAC 3D numerical model with the optimal model parameters to obtain a predicted value of the three-dimensional rockburst crater contour.

[0016] Furthermore, the step S4 of searching for the global optimal parameters includes the following steps:

[0017] Step A1: Set the COA algorithm parameters: number of iterations T, population size N, dimension Dim and temperature Temp; enter the local optimization number I locad ; Maximum number of iterations I max ; Number of crayfish NP;

[0018] Step A2: Set the GPR machine learning algorithm parameters: mean function and covariance function (kernel function) in the activity subset;

[0019] Step A3: Generate an initial random crayfish team NP (i = 1, 2, 3, ..., NP), calculate the optimization objective function value f(X), and sort them from small to large according to the function value to determine the current optimal position X of the food. L ;

[0020] Step A4: Use COA optimization algorithm to perform global optimization, X G represents the optimal position obtained by the number of iterations, X L It represents the optimal position obtained after the previous generation population is updated, that is, the crayfish position and its optimization objective function value, to the historical database X record ;

[0021] Step A5: When the number of local iterations of the COA algorithm is I l The number of times I enters local optimization is reached locad Then, output X L 、f(X L ) and X record ;

[0022] Step A6: Select X record Mid-range X L A certain range of crayfish information is used as training samples;

[0023] Step A7: Use the training samples to train the GPR local proxy model, that is, obtain the original optimization objective function value f(x) in X best The GPR approximate optimization objective function distribution f in the local neighborhood of GPR (x);

[0024] Step A8: Obtain the optimal value f in the distribution of the GPR approximate optimization objective function GPR (x ib ) corresponding to x ib ;

[0025] Step A9: If f GPR (x ib ) is better than f(X L ), then update X L＝ x ib ;

[0026] Step A10: Determine f(X L ) whether the convergence condition is met, if so, then the optimization is terminated; otherwise, return to step A4 until the global iteration number I=I max .

[0027] Furthermore, the fitness function evaluation method in step S2 is to optimize the objective function: the difference between the total number of calculation units within the measured rockburst crater outline and the total number of calculation units that have entered the yield state within the three-dimensional rockburst crater outline calculated by the FLAC 3D numerical simulation model, plus the total number of calculation units that have entered the yield state outside the measured three-dimensional rockburst crater outline multiplied by the penalty factor; its calculation formula is:

[0028] f(X)=N * -N1+N2×C

[0029] Where: N * is the total number of cells divided in the measured damage zone; N1 and N2 are the total number of calculation cells that enter the yield state inside and outside the three-dimensional rockburst crater contour at the end of the FLAC3D numerical simulation; C is the penalty factor, which is 100.

[0030] Furthermore, the training samples of the GPR local proxy model are preferably selected from X record Mid-range X L The location information of the nearest 3 to 5 times the number of crayfish.

[0031] Furthermore, the COA-GPR joint optimization algorithm, through the local optimization times I locad To terminate the COA global optimization, enter the GPR local agent model, select, I max with I locad The ratio should be between 20 and 30.

[0032] Compared with the prior art, the present invention has the following beneficial effects:

[0033] (1) The GPR machine learning algorithm used in the present invention to train the local proxy model has the characteristics of adaptability to small samples, good regression and generalization. Since there is a highly nonlinear mapping relationship between the mechanical parameters of the rock mass and the mechanical response, especially when the dimension of the mechanical parameters in the back analysis is high, using the existing technology, that is, using the machine learning algorithm to train the global proxy model, it is difficult to ensure a certain approximate accuracy requirement under small sample conditions, which can easily cause the problem of optimization failure. On the contrary, in order to improve the proxy accuracy, a large number of training samples are generated, which also brings the problem of excessive numerical reanalysis. Therefore, using the GPR machine learning algorithm to train the local proxy model and perform proxy in the local area of ​​the objective function during the optimization process can reduce the contradiction between numerical reanalysis and proxy accuracy.

[0034] (2) The COA optimization algorithm used in this invention is used for global optimization of the inversion of rock mechanical parameters in deep tunnels. It has the characteristics of high robustness, fast convergence speed, and stronger global optimization ability. Compared with existing algorithms, the COA optimization algorithm does not require complex parameter adjustment or operator design, is easy to implement and apply, and can find more optimal surrounding rock mechanical parameters, making the difference between the measured rockburst crater size parameters and the numerical simulation size parameters smaller, thereby improving the reliability of the rock mechanical parameters.

[0035] (3) The COA-GPR joint optimization algorithm proposed in the present invention is combined with the FLAC3D numerical model for feedback identification of rockburst crater size parameters. The COA optimization algorithm is used to globally optimize the rockburst crater size parameters, and FLAC3D calculates the fitness value. After a certain number of optimization iterations, lobster information or historical lobster information within a certain range of the current optimal temperature in the historical database is selected to establish a training model, and a local proxy model based on GPR is trained to obtain the GPR approximate fitness function distribution of the original fitness function in the local neighborhood, thereby finding a position with better fitness than the current optimal lobster, that is, the difference between the three-dimensional rockburst crater contour currently measured and the three-dimensional rockburst crater contour obtained by FLAC3D numerical simulation is smaller, so that the mechanical parameters of the rockburst crater rock mass that meet the convergence conditions can be quickly found. Compared with the existing single optimization algorithm or the use of a global proxy model, this method effectively improves the efficiency of optimization back analysis while ensuring calculation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 A schematic diagram of the location of rock burst pits in the forward excavation section of a tunnel project provided by the present invention;

[0037] Figure 2 A schematic diagram of the outline of a three-dimensional rock burst pit in a tunnel engineering project provided by the present invention;

[0038] Figure 3 Flowchart based on COA-GPR joint optimization algorithm provided in the present invention

[0039] Figure 4 A flow chart of a collaborative optimization inversion method for mechanical parameters of surrounding rock of a deep tunnel based on a three-dimensional rockburst crater contour is provided in the present invention;

[0040] Figure 5 A schematic diagram of a rockburst simulation result of a tunnel project provided in Example 2 of the present invention;

[0041] Figure 6 Schematic diagram of the food-seeking phase of the crayfish optimization algorithm COA provided in Example 2 of the present invention;

[0042] Figure 7A two-dimensional model diagram of the GPR local optimization experience dataset provided in Example 2 of the present invention;

[0043] Figure 8 The present invention provides a comparison of the number of FLAC3D calls and the search time of two algorithms for searching the global optimal parameters of a deep tunnel project. Specific implementation methods

[0044] The following further illustrates the specific embodiments of the present invention with reference to the accompanying drawings and examples. It should be noted that the accompanying drawings only illustrate portions relevant to the present invention, not all of the embodiments. Furthermore, the specific examples are intended only to illustrate the present invention and are not intended to limit the scope of the invention.

[0045] The present invention provides a method for inverting the mechanical parameters of surrounding rock of a deep tunnel based on a three-dimensional rockburst crater contour and a COA-GPR collaborative optimization algorithm, comprising the following steps:

[0046] Step S1: Constructing a FLAC 3D numerical model of the tunnel: A FLAC 3D numerical simulation model of the tunnel in the area where the rockburst occurs is established. The FLAC 3D numerical simulation model adopts the Mogi-Coulomb hard rock slab crack constitutive model and records the total number and centroid coordinates of the calculation units within the three-dimensional contour of the measured rockburst;

[0047] Step S2: Establish the optimization objective function: Take the mechanical parameters of the tunnel surrounding rock, i.e. the mechanical parameters of the Mogi-Coulomb hard rock slab crack constitutive model, as the optimization variables

[0048] Step S3: Subtract the total number of units in the measured damage zone from the total number of units that have entered the yield state within the three-dimensional rockburst crater contour calculated by FLAC 3D, and add the total number of units that have entered the yield state outside the measured three-dimensional rockburst crater contour multiplied by the penalty factor as the optimization objective function;

[0049] Step S4: Searching for global optimal parameters: Using the crayfish-Gaussian process regression (COA-GPR) collaborative optimization algorithm, combined with the tunnel FLAC 3D numerical model, searching for the optimal solution when the optimization objective function is globally minimized, until the preset value is consistent with the measured value. At this time, the rock mass mechanical parameters in the tunnel surrounding rock mechanical parameter group are the rock mechanical parameter inversion values;

[0050] Step S5: performing forward calculations based on the FLAC 3D numerical model with the optimal model parameters to obtain a predicted value of the three-dimensional rockburst crater contour.

[0051] Furthermore, the step S4 of searching for the global optimal parameters includes the following steps:

[0052] Step A1: Set the COA algorithm parameters: number of iterations T, population size N, dimension Dim and temperature Temp; enter the local optimization number I locad ; Maximum number of iterations I max ; Number of crayfish NP;

[0053] Step A2: Set the GPR machine learning algorithm parameters: mean function and covariance function (kernel function) in the activity subset;

[0054] Step A3: Generate an initial random crayfish team NP (i = 1, 2, 3, ..., NP), calculate the optimization objective function value f(X), and sort them from small to large according to the function value to determine the current optimal position X of the food. L ;

[0055] Step A4: Use COA optimization algorithm to perform global optimization, X G represents the optimal position obtained by the number of iterations, X L It represents the optimal position obtained after the previous generation population is updated, that is, the crayfish position and its optimization objective function value, to the historical database X record ;

[0056] Step A5: When the number of local iterations of the COA algorithm is I l The number of times I enters local optimization is reached locad Then, output X L 、f(X L ) and X record ;

[0057] Step A6: Select X record Mid-range X L A certain range of crayfish information is used as training samples;

[0058] Step A7: Use the training samples to train the GPR local proxy model, that is, obtain the original optimization objective function value f(x) in X best The GPR approximate optimization objective function distribution f in the local neighborhood of GPR (x);

[0059] Step A8: Obtain the optimal value f in the distribution of the GPR approximate optimization objective function GPR (x ib ) corresponding to x ib ;

[0060] Step A9: If f GPR (x ib ) is better than f(X L ), then update X L =x ib ;

[0061] Step A10: Determine f(X L ) whether the convergence condition is met, if so, then the optimization is terminated; otherwise, return to step A4 until the global iteration number I=I max .

[0062] Furthermore, the fitness function evaluation method in step 2 is to optimize the objective function: the difference between the total number of calculation units within the measured rockburst crater outline and the total number of calculation units that enter the yield state within the three-dimensional rockburst crater outline calculated by the FLAC 3D numerical simulation model, plus the total number of calculation units that enter the yield state outside the measured three-dimensional rockburst crater outline multiplied by the penalty factor; its calculation formula is:

[0063] f(X)=N * -N1+N2×C

[0064] Where: N * is the total number of cells divided in the measured damage zone; N1 and N2 are the total number of calculation cells that enter the yield state inside and outside the three-dimensional rockburst crater contour at the end of the FLAC3D numerical simulation; C is the penalty factor, which is 100.

[0065] Furthermore, the training samples of the GPR local proxy model are preferably selected from X record Mid-range X L The location information of the nearest 3 to 5 times the number of crayfish.

[0066] Furthermore, the COA-GPR joint optimization algorithm, through the local optimization times I locad To terminate the COA global optimization, enter the GPR local agent model, select, I max with I locad The ratio should be between 20 and 30.

[0067] The Mogi-Coulomb constitutive model for hard rock slab cracking in step 1 is described as follows:

[0068] (1) Mogi-Coulomb model

[0069] The Mogi-Coulomb criterion considers the octahedral strength of the intermediate principal stress and can effectively characterize the nonlinear mechanical behavior of deeply buried rock masses under high stress, offering significant advantages in describing the mechanical behavior of deeply buried, high-stress rock masses. Furthermore, the Mogi-Coulomb criterion addresses the issue of discontinuities in the derivatives of the yield surface of the traditional Mohr-Coulomb criterion at the corners of a hexagonal cone in space, in principal stress space.

[0070] τ oct =a+bσ m,2 (1)

[0071]

[0072]

[0073] Where: τ oct is the octahedral shear stress, σ m,2 is the average effective stress acting on the shear surface, a and b are parameters related to rock medium properties, corresponding to the intercept and slope, respectively.

[0074] When the intermediate principal stress is equal to the third principal stress, that is, σ2 = σ3, the Mogi-Coulomb criterion degenerates into the Mohr-Coulomb criterion, and we can get:

[0075]

[0076] Where: c, are the cohesion and internal friction angle corresponding to the rock medium, respectively.

[0077] Because the Mogi-Coulomb constitutive model's description of the post-peak mechanical behavior of deep-buried surrounding rock is significantly different from the actual mechanical behavior without considering elastic-plastic damage, elastic-plastic damage correction is required to reasonably describe the nonlinear mechanical behavior of deep-buried tunnel surrounding rock. The Mogi-Coulomb yield surface function and plastic potential surface function are expressed as:

[0078]

[0079]

[0080]

[0081] Where: J2 is the second invariant of the deviatoric stress tensor; S is the damage variable; a', b' are rock material parameters, φ is the expansion angle.

[0082] The evolution process of rock damage variables is characterized by equivalent plastic strain, and the equivalent plastic strain ε is introduced. q To characterize the evolution of rock damage variables, the damage variable S is an exponential function of the equivalent plastic strain.

[0083]

[0084] Where: ε 1p , ε 2p and ε 3p There are three principal plastic strains.

[0085] The rock compression process is divided into the undamaged and damaged stages, and the corresponding evolution equation of the damage variable S is:

[0086]

[0087] Where: η is the parameter obtained from the experimental fitting, reflecting the evolution rate of the damage variable S with plastic strain; is the equivalent strain when damage occurs.

[0088] when When , that is, when the rock internal friction angle is equal to the rock expansion angle, the flow law associated with the plastic potential surface can be obtained as follows:

[0089]

[0090] Where: is the plastic strain increment, dλ is the plasticity factor and dλ≥0.

[0091] The non-correlated flow law is:

[0092]

[0093] Where: N is the flow vector related to the plastic potential surface:

[0094]

[0095] Where: C1, C2 are the flow vector parameter coefficients related to the plastic potential surface, Substituting into formula (12) we get:

[0096]

[0097] Where: S x 、S y 、S z , τ yz , τ zx , τ xy are the stress components.

[0098] At the same time, considering the material damage:

[0099]

[0100] Where: is the elastic strain tensor; is the elastic stiffness matrix of the material after damage.

[0101] Under the isotropic condition of rock medium, by reducing the elastic stiffness matrix, the damage effect of the material and the plastic flow can be reflected to describe the strain softening stage of the material, namely:

[0102]

[0103] Where: De is the elastic stiffness matrix.

[0104] where ε ps is the irrecoverable strain caused by material damage and plastic flow, which includes damage strain, ε es It is the elastic strain generated after the material is damaged and yielded.

[0105] For the modified constitutive model, if it is to be developed and applied in numerical simulation, since the numerical simulation software mostly calculates the constitutive model in an incremental iterative form, it is necessary to establish an incremental constitutive relationship before the constitutive model can be applied to numerical simulation calculations. Therefore, the incremental constitutive relationship based on Mogi-Coulomb elastic-plastic damage is further derived:

[0106] (1) When the rock is in the elastic stage, the constitutive relation obeys the generalized Hooke's law:

[0107]

[0108] (2) When the rock is in the plastic stage, the total strain increment of the rock is the sum of the elastic strain increment and the plastic strain increment. Based on this, the elastic-plastic damage constitutive model of the rock is:

[0109]

[0110] Substituting Equation (10) into Equation (17), the non-associated Mogi-Coulomb elastic-plastic damage constitutive equation can be obtained as follows:

[0111]

[0112] When the material is an ideal elastic-plastic material, the consistency condition is satisfied when plastic deformation occurs, that is, the stress point is always on the yield surface during plastic deformation, and the following conditions are satisfied:

[0113]

[0114] Combining Equations (11) and (19) and substituting them into Equation (20) allows us to solve the plasticity factor λd. This completes the establishment of the Mogi-Coulomb elastic-plastic damage increment constitutive relation. The above principle method is encapsulated according to the requirements of the FLAC3D constitutive model and used for FLAC3D numerical calculations.

[0115] It should be noted that establishing a new constitutive model in FLAC 3D follows the same operational principles as the original constitutive model: given the stress state at time t and the total strain increment corresponding to the Δt time step, the stress increment corresponding to the Δt time step and the updated stress state at time t+Δt are calculated. Therefore, the Mogi-Coulomb model and constitutive model form for hard and brittle tunnel rock under high geostress are adopted.

[0116] The tunnel surrounding rock mechanical parameters described in step S2 are used as optimization variables to construct a minimized optimization objective function, which is specifically described as follows:

[0117] The optimal parameter combination is searched through iterative calculation, so that the total number of units in the measured damage zone and the total number of units that have entered the yield state within the three-dimensional rockburst crater contour calculated by FLAC 3D are subtracted, and the total number of units that have entered the yield state outside the measured three-dimensional rockburst crater contour multiplied by the penalty factor is used as the optimization objective function. The parameters to be determined are gradually approximated and corrected. The objective function can be set as:

[0118] f(X)=N * -N1+N2×C (20)

[0119] Where: N * is the total number of cells divided in the measured damage zone; N1 and N2 are the total number of calculation cells that enter the yield state inside and outside the three-dimensional rockburst crater contour at the end of FLAC3D numerical simulation; C is the penalty factor, which is generally taken as 100.

[0120] The specific description of the COA-GPR joint optimization algorithm in step S4 is as follows:

[0121] (1) The crayfish optimization algorithm is an emerging meta-heuristic optimization algorithm. The basic implementation steps of the algorithm are mainly divided into: crayfish foraging, summer escape and competition behavior. Its principle is as follows:

[0122] ①Algorithm initialization

[0123] Initialize a group of crayfish individuals, each of which has a position in the search space, which represents a possible solution to the optimization problem. These individuals have different properties, such as position, speed, and sensing range.

[0124] ② Location update

[0125] The crayfish's position update is usually based on the current position, the best position, and a certain amount of random perturbation. In each iteration, the crayfish individuals update their positions according to the following factors:

[0126] a. Fitness of the current location: This measures the quality of the current solution. If the fitness of the current location is good, the individual will tend to continue exploring near that area.

[0127] b. Interactions between individuals: Each crayfish is influenced by other individuals within a certain range. If there is a more optimal individual nearby, it will move a certain distance towards that individual.

[0128] c. Random perturbation: In order to maintain diversity, the algorithm will also introduce random perturbations.

[0129] ③ Stage division

[0130] a. Summer Escape (Exploration): When temperatures exceed 30°C, crayfish will choose to escape the heat in burrows. In the algorithm, this represents the process of searching for a superior solution. If there are no other crayfish competing for the burrow, the crayfish will directly enter the burrow; otherwise, competition will occur.

[0131] b. Competition stage (development stage): When there are other crayfish competing for the same burrow, they will compete with each other and fight for the burrow by adjusting their positions.

[0132] c. Foraging (Development): When the temperature is suitable, crayfish begin foraging. In the algorithm, this represents the process of further developing optimal solutions. Depending on the size of the food, crayfish will choose to eat it directly or tear it into pieces before eating.

[0133] In summary, the optimization process of the COA algorithm is as follows:

[0134] Step (1), initialize the population, calculate the fitness value of the population and obtain X G and X L .

[0135] Step (2) defines the living environment of crayfish according to equations (21) and (22).

[0136] Step (3): When the temperature is greater than 30 degrees and rand < 0.5, COA obtains a new position according to equation (26) and enters step 8.

[0137] Step (4): When the temperature is greater than 30 degrees and rand ≥ 0.5, COA obtains a new position according to equations (28) and (29) and enters step 8.

[0138] Step (5): When the temperature is less than or equal to 30°C, COA enters the foraging phase and defines the food intake p and food size Q according to equations (30) and (31).

[0139] Step (6): If Q > (C3 + 1) / 2, shred the food according to equation (32). Then, obtain a new position by ingesting the food according to equation (33) and proceed to step 8.

[0140] Step (7), if Q≤(C3+1) / 2, obtain the new position through equation (34) and go to step 8.

[0141] Step (8) evaluates whether the population has exited the loop. If not, return to step 2.

[0142] Step (9): output the individual with the best position.

[0143] Preferably, in the initialization phase of the multidimensional optimization problem, each crayfish represents a 1×dim matrix, with each column of the matrix representing a solution to the problem. COA is initialized by randomly generating N sets of candidate solutions X between upper and lower bounds. N is the population size, and dim is the population dimension. COA is initialized as follows:

[0144] X=[X1,X2,....,X N ] (twenty one)

[0145] X ij =lb j +(ub j -lb j )*rand (22)

[0146] X ij is the position of individual i in dimension j, where lb j represents the lower bound of the j-th dimension, ub j represents the upper bound of the j-th dimension, and rand is a random number in [0,1].

[0147] Step S4.1: Define temperature and crayfish feeding amount

[0148] Changes in temperature will affect the behavior of crayfish, causing them to enter different stages. The definition of temperature is shown in equation (3). When the temperature exceeds 30°C, crayfish will choose a cool place to escape the heat. At the appropriate temperature, crayfish will engage in foraging behavior. The amount of food consumed by crayfish is affected by temperature. The feeding range of crayfish is between 15 and 30°C, with 25°C being the best. Therefore, the amount of food consumed by crayfish can be approximated to a normal distribution, so that the amount of food consumed is affected by temperature. The mathematical model of crayfish food intake and the amount of food consumed at different temperatures are as follows Figure 6 As shown:

[0149] temp=rand×15+20 (23)

[0150] Here, temp represents the temperature of the environment where the crayfish is located.

[0151]

[0152] Among them, μ refers to the temperature that is most suitable for crayfish, and σ and C1 are used to control the intake of crayfish at different temperatures.

[0153] Step S4.2: Summer vacation phase (exploration phase)

[0154] When the temperature is above 30 degrees, it means the temperature is too high. At this time, crayfish will enter the cave to escape the heat. The definition of the cave is as follows:

[0155] X shade =(X G +X L ) / 2 (25)

[0156] where X G represents the optimal position obtained by the number of iterations, X L Indicates the optimal position obtained after the previous generation population update.

[0157] The competition among crayfish for the burrow is a random event. In COA, when rand < 0.5, it means that there are no other crayfish competing for the burrow, and the crayfish will directly enter the burrow to avoid the heat. The formula for crayfish entering the burrow to avoid the heat is as follows:

[0158]

[0159] Where t represents the current iteration number, t+1 represents the next iteration number, and C2 is the decreasing curve.

[0160] C2=2-(t / T) (27)

[0161] Where T represents the maximum number of iterations.

[0162] Step S4.3: Competition phase (development phase)

[0163] When the temperature is greater than 30 degrees Celsius, rand ≥ 0.5. This means that other crayfish have also chosen this burrow. At this point, they compete for the burrow. They compete for the burrow using the following formula.

[0164]

[0165] z=round(rand*(N-1))+1 (29)

[0166] Where z represents a random individual crayfish.

[0167] Step S4.4: Foraging phase (development phase)

[0168] When the temperature is less than or equal to 30 degrees, the temperature is suitable for crayfish to eat. At this time, crayfish will go out to look for food. When eating, crayfish will choose whether to tear the food into pieces according to the size of the food. If the food is of the right size, the crayfish will use the claws to tear the food into pieces and then use the second and third walking legs to alternately pick up the food and eat it. Food is defined as:

[0169] X food =X G (30)

[0170] Food size is defined as:

[0171] Q=C3*rand*(fitness i / fitness food ) (31)

[0172] C3 is the food factor, which represents the maximum food and has a value of constant 3. i Represents the fitness value of the i-th crayfish, fitness food Represents the fitness value of the location where the food is located.

[0173] When Q>(C3+1) / 2, it means the food is too big. At this time, the crayfish will tear the food into pieces according to the following formula.

[0174]

[0175] After tearing the food into pieces, the crayfish will alternately use its second and third walking legs to pick up the food and ingest it. To simulate this alternating feeding behavior, a combination of sine and cosine functions is used in the equation to simulate the alternating process, as shown in the figure. Furthermore, the amount of food a crayfish obtains is also related to the amount of food it ingests. The feeding equation is as follows:

[0176]

[0177] When Q≤(C3+1) / 2, the crayfish will move directly toward the food and eat it. The equation is as follows:

[0178]

[0179] (2) GPR is a machine learning model suitable for solving complex regression problems with small samples, high dimensions, and nonlinearity. Compared with machine learning models such as SVM and ANN, GPR is easier to implement and has significant advantages such as maximum probability prediction output and parameter adaptation. All statistical characteristics of GPR are determined by the Gaussian distribution of its mean m(t) and covariance function k(t,t'):

[0180] f(t)GP(m(t),k(t,t')) (35)

[0181] Assume that a training sample set has n observation data, and the training sample set can be expressed as D = {(x i ,y i |x i =1,…,n)}, x is the d-dimensional input vector, and the observed target value y i ∈R. If X represents a d×n dimensional input matrix and y represents an output vector, then the training sample set can be expressed as D=(X,y). The GPR model is trained using the training samples, and thus the output vector corresponding to the new input X can be predicted.* The corresponding output y * .

[0182] Assuming that X follows a Gaussian distribution, the linear regression model is:

[0183] y=f(X)+ε (36)

[0184] Where f(X) is the regression function value, the noise ε is in accordance with the Gaussian distribution, with a mean of 0 and a variance of σ n 2, namely:

[0185]

[0186] The prior distribution of the observed target value y is:

[0187]

[0188] Where: K = K(X,X) is an n×n symmetric positive definite covariance matrix, K ij Representation x i with x j correlation.

[0189] The joint Gaussian prior distribution of the regression function output of the training sample and the test sample is:

[0190]

[0191]

[0192] Among them, K(X,X * ) represents the test input point X * With the n×1 order covariance matrix of all input points X in the training set; k(X * ,X * ) represents the test input point X * The covariance function used is:

[0193]

[0194] Among them, l, σ f , σ n is a hyperparameter that can be obtained by the maximum likelihood method.

[0195] According to the Bayesian principle, the posterior probability is deduced and X * The corresponding output y * The predicted mean and variance of :

[0196]

[0197]

[0198] Once the distribution of f(x) is determined, the function value of f(x) at x can be predicted, which is the prediction process.

[0199] (3) The flowchart of COA-GPR joint optimization algorithm is shown in Figure 3 , the specific implementation steps are as follows:

[0200] Step A1: Set the COA algorithm parameters: number of iterations T, population size N, dimension Dim and temperature Temp; enter the local optimization number I locad ; Maximum number of iterations I max ; Number of crayfish NP;

[0201] Step A2: Set the GPR machine learning algorithm parameters: mean function and covariance function (kernel function) in the activity subset;

[0202] Step A3: Generate an initial random crayfish team NP (i = 1, 2, 3, ..., NP), calculate the optimization objective function value f(X), and sort them from small to large according to the function value to determine the current optimal position X of the food. L ;

[0203] Step A4: Use COA optimization algorithm to perform global optimization, X G represents the optimal position obtained by the number of iterations, X L It represents the optimal position obtained after the previous generation population is updated, that is, the crayfish position and its optimization objective function value, to the historical database X record ;

[0204] Step A5: When the number of local iterations of the COA algorithm is I l The number of times I enters local optimization is reached locad Then, output X L 、f(X L ) and X record ;

[0205] Step A6: Select X record Mid-range X L A certain range of crayfish information is used as training samples;

[0206] Step A7: Use the training samples to train the GPR local proxy model, that is, obtain the original optimization objective function value f(x) in X best The GPR approximate optimization objective function distribution f in the local neighborhood of GPR (x);

[0207] Step A8: Obtain the optimal value f in the distribution of the GPR approximate optimization objective function GPR (x ib ) corresponding to x ib ;

[0208] Step A9: If f GPR (x ib ) is better than f(X L ), then update X L =x ib ;

[0209] Step A10: Determine f(X L ) whether the convergence condition is met, if so, then the optimization is terminated; otherwise, return to step A4 until the global iteration number I=I max .

[0210] (4) The COA-GPR joint optimization algorithm is used, combined with the feedback identification steps of the 3D rockburst pit FLAC3D numerical model as follows:

[0211] ① Select n model parameters from the numerical model parameters as optimization variables, determine their optimization range, and select other model parameters from the original value range; determine the convergence condition of the optimization objective function f(X) to be less than the threshold f Best ;

[0212] ② Based on the COA-GPR joint optimization algorithm, determine the COA algorithm parameters: number of iterations T, population size N, dimension Dim and temperature Temp; enter the local optimization number I locad ; Maximum number of iterations I max ; Number of crayfish NP;

[0213] ③ Generate an initial random team NP (i = 1, 2, 3, ..., NP), and perform FLAC3D forward calculation on the crayfish group to obtain the optimization objective function value f(X), and sort them from small to large according to the fitness function value to determine the current optimal position X of the food. L ;

[0214] ④ Use the COA optimization algorithm to perform global optimization and record the crayfish position information of each iteration, that is, the crayfish position and the objective function value calculated by FLAC3D, to the historical database X record ;

[0215] ⑤ When the number of local iterations of the COA algorithm is I l The number of times I enters local optimization is reached locad Then, output X L 、f(X L ) and X record ;

[0216] ⑥Select X record Mid-range X L The location information of the nearest 3 to 5 times the number of crayfish is used as training samples;

[0217] ⑦ Use the training samples to train the GPR local proxy model, that is, to obtain the original optimization objective function value f(x) in X best The GPR approximate optimization objective function distribution f in the local neighborhood of GPR (x);

[0218] ⑧ Obtain the objective function value through FLAC3D forward calculation, if f GPR (x ib ) is better than f(X L ), then update X L＝ x ib ;

[0219] ⑨If f(X L ) <f Best , then the optimization is ended and the optimal numerical model parameters are output; otherwise, return to step ④ until the number of iterations I=I max .

[0220] The COA-GPR joint optimization algorithm is selected through the local optimization number I locad To terminate the COA global optimization and enter the GPR local proxy model, I max with I locad The ratio should be between 20 and 30.

[0221] Example 1

[0222] To understand the performance of the COA-GPR joint optimization algorithm, Example 1 compares the performance of the COA-GPR joint optimization algorithm proposed in this paper with the particle swarm optimization algorithm (PSO), the gray wolf optimization algorithm (GWO), and the moth flame optimization algorithm (MFO) based on some CEC2005 benchmark test functions. The algorithm parameter settings are shown in Table 1. The superiority of the algorithm is verified by the optimal value, average value, and standard deviation.

[0223] Table 1 Algorithm parameter settings

[0224]

[0225] CEC2005 benchmark test function, where f1-f7 are high-dimensional unimodal functions, f8-f 13 is a high-dimensional multimodal function, f 14 —f 23 is a fixed-dimensional multimodal function. The high-dimensional unimodal functions of f1-f7 refer to functions with only one minimum value. The functions listed here can test the global development performance of the algorithm. 13 A high-dimensional multimodal function is one that has multiple extreme values ​​and increases with the increase of dimension. This function can test the local exploration performance of the algorithm. 14 —f 23It is a fixed-dimensional multimodal function, and its most obvious feature is that the dimension is fixed. We select f2, f6, f9, and f 11 , f 16 , f 21 , see Table 2.

[0226] As can be seen from Table 2, the convergence accuracy of the COA-GPR algorithm is better than that of the PSO algorithm, GWO algorithm, and MFO algorithm, and the global optimization ability is stronger. This shows that the COA-GPR joint optimization algorithm has significant advantages in typical test functions.

[0227] Table 2 Optimization results based on different algorithms

[0228]

[0229] Example 2

[0230] See also Figure 5 、 Figure 6 、 Figure 7 and Figure 8 As shown in the figure, during the bottom construction of the 9+197 to 9+212 section of a deep-buried water diversion tunnel (constructed by the bench method), a strong rockburst occurred at the junction of the upper and lower steps of the tunnel, resulting in the failure of the support in this section, the collapse of the side wall rock mass, and the formation of a large-scale rockburst cavity. The specific implementation steps are as follows:

[0231] The geometric contours of the rockburst pit were measured on site, and the inversion method of the mechanical parameters of the tunnel surrounding rock was carried out on this basis. The specific implementation process is shown in Figure 3 A new intelligent algorithm based on the Crayfish Optimization Algorithm (COA) and Gaussian Process Regression (GPR) machine learning combines the low-computational cost of the GPR proxy model with the powerful global optimization capabilities of the COA. This approach, combined with the numerical calculation software FLAC3D, proposes the COA-GPR-FLAC3D method for inversion of rock mass mechanical parameters. The specific implementation steps are as follows:

[0232] ① Initialize the algorithm parameters: COA algorithm parameters: N = 80, dim = 2.0, temp = 25°, dimension upper and lower bounds lb j and ub j =[1,1]; the initial hyperparameters of the GPR model are: lni=[-1,-1,-1,-1,-1,-1,-1], lnσ f =-1, lnσ n =ln1×10 -6After 7 iterations of the COA algorithm, GPR machine learning is started. The total number of learning samples selected in the neighborhood of the local optimal individual is 3×80=240. GPR machine learning selects 100 of the 240 learning samples as information vectors each time, and repeats the selection process 8 times to obtain the 100 most "useful" information vectors for learning regression. The convergence criterion of the objective function of the three algorithms is ε=1×10 -3 ;

[0233] ② Set the range of the algorithm’s individual search space, which is the range of values ​​of the surrounding rock mechanical parameters to be inverted (see Table 3);

[0234] Table 3 Search intervals for surrounding rock mechanical parameters

[0235]

[0236] ③ Start the optimization algorithm to perform optimization operation. After the individual evolves according to the algorithm's evolution strategy, the new position of the individual (80 sets of surrounding rock mechanical parameters) is saved in the data interface file A;

[0237] ④ Start the FLAC 3D numerical calculation software through the custom software call command, read the surrounding rock mechanical parameters in the interface file A through the FISH program embedded in FLAC 3D, and substitute them into the established FLAC 3D numerical model to obtain the calculated N*, which is the total number of units divided in the measured damage zone; N1 and N2 are the three-dimensional rock burst pit contours at the end of the FLAC3D numerical simulation (see Figure 1 and Figure 2 The total number of calculation units that enter the yield state inside and outside is shown, and the objective function f(X)=N * -N1+N2×C to obtain the objective function values ​​corresponding to 80 sets of surrounding rock mechanical parameters. The penalty factor C is used in GPR to control the degree of penalty for out-of-cell classification errors in 3D contours to improve the generalization ability of the model. The objective function values ​​are stored in data file B.

[0238] ⑤ Use MATLAB program to read the objective function value in data file B, obtain the fitness value of all individuals, and select the minimum fitness value by comparison. The position coordinates of the individual corresponding to the minimum fitness value are the optimal surrounding rock parameter combination corresponding to this generation of individuals;

[0239] ⑥ Combine this minimum fitness value with the convergence criterion ε=1×10 -3 Compare and if it is less than the convergence criterion, stop the calculation and output the individual position coordinates (i.e. the optimal combination of surrounding rock mechanical parameters); otherwise, continue a new round of optimization calculation until the convergence criterion is reached;

[0240] The parameter inversion results of the two optimization algorithms are shown in Table 4. Compared with COA, the parameter inversion results of the COA-GPR algorithm are closer to the true values;

[0241] Table 4 Inversion results of surrounding rock mechanical parameters

[0242]

[0243] The cost and computational time of parameter inversion using different algorithms are shown in Table 5. It can be seen that the computational cost and computational time of COA-GPR are much lower than those of COA, indicating that the COA-GPR collaborative optimization algorithm has significant advantages in the inversion of mechanical parameters of surrounding rock masses of large cavern groups.

[0244] Table 5 Comparison of computational time of two algorithms

[0245] algorithm COA COA-GPR FLAC3D call count 2660 1504 Calculation time (s) <![CDATA[4.18×10 3 ]]> <![CDATA[3.12×10 3 ]]>

[0246] It should be noted that the purpose of disclosing the above examples is to facilitate a further understanding of the present invention. However, those skilled in the art will appreciate that various obvious changes, readjustments, and substitutions to the present invention do not depart from the scope of protection of the present invention. Therefore, the present invention is not limited to the contents disclosed in the examples, and the scope of protection claimed by the present invention shall be determined by the scope defined in the claims.

Claims

1. A method for inversion of mechanical parameters of surrounding rock of deep buried tunnels based on three-dimensional rockburst crater contour and COA-GPR collaborative optimization algorithm, characterized by: The steps include: Step S1: Constructing a FLAC 3D numerical model of the tunnel: A FLAC 3D numerical simulation model of the tunnel in the area where the rockburst occurs is established. The FLAC 3D numerical simulation model adopts the Mogi-Coulomb hard rock slab crack constitutive model and records the total number and centroid coordinates of the calculation units within the three-dimensional contour of the measured rockburst; Step S2: Establish the optimization objective function: Take the mechanical parameters of the tunnel surrounding rock, i.e. the mechanical parameters of the Mogi-Coulomb hard rock slab crack constitutive model, as the optimization variables Step S3: Subtract the total number of units in the measured damage zone from the total number of units that have entered the yield state within the three-dimensional rockburst crater contour calculated by FLAC 3D, and add the total number of units that have entered the yield state outside the measured three-dimensional rockburst crater contour multiplied by the penalty factor as the optimization objective function; Step S4: Searching for global optimal parameters: Using the crayfish-Gaussian process regression (COA-GPR) collaborative optimization algorithm, combined with the tunnel FLAC 3D numerical model, searching for the optimal solution when the optimization objective function is globally minimized, until the preset value is consistent with the measured value. At this time, the rock mass mechanical parameters in the tunnel surrounding rock mechanical parameter group are the rock mechanical parameter inversion values; Step S5: performing forward calculations based on the FLAC 3D numerical model with the optimal model parameters to obtain a predicted value of the three-dimensional rockburst crater contour.

2. The method for inversion of mechanical parameters of surrounding rock of deep tunnels based on three-dimensional rockburst crater contour and COA-GPR collaborative optimization algorithm according to claim 1 is characterized in that: The step S4 of searching for the global optimal parameters comprises the following steps: Step A1: Set the COA algorithm parameters: number of iterations T, population size N, dimension Dim and temperature Temp; enter the local optimization number I locad ; Maximum number of iterations I max ; Number of crayfish NP; Step A2: Set the GPR machine learning algorithm parameters: mean function and covariance function (kernel function) in the activity subset; Step A3: Generate an initial random crayfish team NP (i = 1, 2, 3, ..., NP), calculate the optimization objective function value f(X), and sort them from small to large according to the function value to determine the current optimal position X of the food. L ; Step A4: Use COA optimization algorithm to perform global optimization, X G represents the optimal position obtained by the number of iterations, X L It represents the optimal position obtained after the previous generation population is updated, that is, the crayfish position and its optimization objective function value, to the historical database X record ; Step A5: When the number of local iterations of the COA algorithm is I l The number of times I enters local optimization is reached locad Then, output X L 、f(X L ) and X record ; Step A6: Select X record Mid-range X L A certain range of crayfish information is used as training samples; Step A7: Use the training samples to train the GPR local proxy model, that is, obtain the original optimization objective function value f(x) in X best The GPR approximate optimization objective function distribution f in the local neighborhood of GPR (x); Step A8: Obtain the optimal value f in the distribution of the GPR approximate optimization objective function GPR (x ib ) corresponding to x ib ; Step A9: If f GPR (x ib ) is better than f(X L ), then update X L =x ib ; Step A10: Determine f(X L ) whether the convergence condition is met, if so, then the optimization is terminated; otherwise, return to step A4 until the global iteration number I=I max .

3. The method for inversion of mechanical parameters of surrounding rock of deep tunnels based on three-dimensional rockburst crater contour and COA-GPR collaborative optimization algorithm according to claim 1 is characterized in that: The fitness function evaluation method in step S2 is to optimize the objective function: the difference between the total number of calculation units within the measured rockburst crater outline and the total number of calculation units that have entered the yield state within the three-dimensional rockburst crater outline calculated by the FLAC 3D numerical simulation model, plus the total number of calculation units that have entered the yield state outside the measured three-dimensional rockburst crater outline multiplied by the penalty factor; its calculation formula is: f(X)=N * -N1+N2×C Where: N * is the total number of cells divided in the measured damage zone; N1 and N2 are the total number of calculation cells that enter the yield state inside and outside the three-dimensional rockburst crater contour at the end of the FLAC3D numerical simulation; C is the penalty factor, which is 100.

4. The method for inversion of mechanical parameters of surrounding rock of deep tunnels based on three-dimensional rockburst crater contour and COA-GPR collaborative optimization algorithm according to claim 2, characterized in that: The training samples of the GPR local proxy model are preferably selected as X record Mid-range X L The location information of the nearest 3 to 5 times the number of crayfish.

5. The method for inversion of mechanical parameters of surrounding rock of deep tunnels based on three-dimensional rockburst crater contour and COA-GPR collaborative optimization algorithm according to claim 1, characterized in that: COA-GPR joint optimization algorithm, through the local optimization times I locad To terminate the COA global optimization, enter the GPR local agent model, select, I max with I locad The ratio should be between 20 and 30.

Images