Deep tunnel surrounding rock mechanical parameter inversion method based on three-dimensional brittle failure zone contour and SSA-IVM joint optimization algorithm

Abstract

The invention discloses a deep tunnel surrounding rock mechanical parameter inversion method based on a three-dimensional brittle failure zone contour and sparrow optimization algorithm-information vector machine (SSA-IVM) combined optimization algorithm. The engineering technical problem that due to the fact that excavation instantaneous displacement is difficult to monitor, deep tunnel surrounding rock mechanical parameters are difficult to reasonably obtain through displacement back analysis is solved. The method comprises the following steps: firstly, constructing a tunnel FLAC3D numerical simulation model with the same ground stress condition at the occurrence position of a three-dimensional brittle failure zone; secondly, taking an absolute error between the total number of computational grid units in the actually measured brittle failure zone and the total number of computational grid units entering a plastic state in the range of the actually measured brittle failure zone after calculation of the FLAC3D numerical model as an optimization objective function; and then, by taking the tunnel surrounding rock mechanical parameters as optimization variables and taking a target function reaching a global minimum value as a target, performing global optimization by combining a tunnel FLAC3D numerical model and adopting an SSA-IVM joint optimization algorithm, thereby obtaining reasonable tunnel surrounding rock mechanical parameters.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of rock engineering, and particularly relates to a deep-buried tunnel surrounding rock mechanical parameter inversion method based on a three-dimensional brittle failure zone profile and a SSA-IVM combined optimization algorithm. BACKGROUND

[0002] The rigid demand of China in economic construction, social development, resource utilization, national security and other aspects is increasing, and the underground space gradually becomes a new target for our exploration, utilization and development. However, in the process of the continuous development of underground engineering construction to the deep part, the brittle failure disaster is more prominent, which seriously endangers the production and life safety.

[0003] Brittle failure, as a typical form of geological disasters in underground engineering, is more likely to occur in geological weak areas such as cleavage planes or structural planes of rock. The geological structure of these areas is complex, making it difficult to directly and accurately determine the mechanical properties of surrounding rock. One of the important bases for this classification is the extension depth of the brittle failure zone, which can reflect the scale of brittle failure to a certain extent. However, compared with the depth, the damaged surface area and the volume of the failure zone caused by brittle static failure, as a multi-dimensional measurement index, can more comprehensively and reasonably characterize the damage scale. The larger the rock mass damage surface area and the failure zone scale, the more rock blocks are stripped, and the consequences and losses caused by them are often more serious. Therefore, it is necessary to use a new method, that is, an optimization inverse analysis technology based on the measured data of the brittle failure zone, that is, after the brittle static failure occurs, the measured failure zone size is analyzed to obtain accurate numerical model parameters. Then, using these numerical models with accurate parameters to simulate the subsequent tunnel excavation process can provide scientific and strong basis for the reasonable support design and safety stability evaluation of deep-buried tunnel surrounding rock.

[0004] For the existing optimization inverse analysis method, global optimization algorithms such as genetic algorithm, particle swarm optimization algorithm, evolutionary difference algorithm and ant colony algorithm are often used for global optimization, but it needs to call the time-consuming tunnel numerical calculation model for thousands of times or even hundreds of thousands of times, and the inverse analysis efficiency is low. In practical application, the total calculation time is too long, and it is not suitable for optimization inverse analysis of tunnel numerical model parameters. Combining efficient optimization algorithm, machine learning and numerical calculation can significantly reduce the number of calls to the numerical calculation model, and significantly improve the efficiency of inverse analysis.

[0005] Sparrow optimization algorithm is a swarm intelligence optimization algorithm inspired by the foraging behavior and anti-predation behavior of sparrows in nature. It is a new type of optimization algorithm proposed in 2021, which is inspired by the group behavior and strategy exhibited by sparrows during foraging. When searching for food, sparrows use a discoverer-follower model and add a detection and warning mechanism to ensure the safety and efficiency of the foraging process. Sparrow optimization algorithm is an optimization algorithm with high convergence accuracy, fast convergence speed, strong robustness and strong global search ability.

[0006] Informative Vector Machine (IVM) is an approximate algorithm based on Gaussian process algorithm. Due to the use of Bayesian statistical learning theory and kernel method, IVM has many advantages such as adaptive super parameter acquisition, strong adaptability to high dimension and complex nonlinear problems, and probability significance of prediction output. At the same time, IVM maintains the process variance implied by the kernel function, which can be tracked and used for active subset selection to provide a sparse representation for the model, thereby significantly reducing the time and space complexity of learning. Compared with Gaussian process, IVM algorithm can reduce computational complexity and memory usage without losing accuracy.

[0007] The basic idea of the collaborative optimization algorithm based on SSA-IVM is to use the strong optimization ability of SSA algorithm to perform global optimization (foraging) and record the position information of the optimization operator (sparrow). After a certain number of iterations (after searching for food for a period of time), select the operators (other sparrows or their footprints) within a certain range of the optimization process (sparrow footprint) and the current optimal operator (the current sparrow's most suitable foraging path and outside the range of the predator) to form a training sample, train the IVM local proxy model, and build an approximate expression of the target function (food) in the local range of the current optimal operator (the current sparrow's most suitable foraging path and outside the range of the predator). Thus, the IVM local proxy model is used to predict the optimal operator compared with the optimal operator of SSA global optimization, and the optimal operator is updated to achieve convergence (successful foraging). SUMMARY

[0008] The present application aims to, in high stress tunnel excavation, fault, joint, bedding and other geological structures not only change the continuity of rock, but also significantly reduce its mechanical properties and stability, leading to the easy instability and brittleness of surrounding rock under disturbance. The complexity of these geological structure surfaces makes the stress distribution of rock uneven, and the damage mode localized, resulting in the difficulty of accurately reflecting the true mechanical behavior of rock by the macroscopic mechanical parameters. An intelligent joint optimization inversion method for deep buried tunnel surrounding rock mechanical parameters based on three-dimensional brittle failure zone profile is proposed, which provides a scientific basis for the reasonable evaluation of surrounding rock stability and the reasonable support design in the excavation process of deep buried underground tunnel engineering.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] Step 1: Establish the FLAC3D numerical model of the tunnel. A hard rock constitutive model is used to establish the FLAC3D numerical model of the tunnel.

[0011] Step 2: Establish the optimization objective function. The objective function is the degree of agreement between the measured brittle fracture zones and the simulated brittle fracture zones output by FLAC3D.

[0012] Step 3: Search for globally optimal parameters. After brittle fracture occurs, based on the measured values ​​of the brittle failure zone of the tunnel sidewall or face, the SSA-IVM joint optimization algorithm is used in conjunction with the tunnel FLAC3D numerical model to search for the optimal model rock mechanics parameters until the preset value matches the measured value. At this point, the rock mechanics parameters in the rock mechanics parameter set are the inverted values ​​of the rock mechanics parameters.

[0013] Step 1 detailed explanation:

[0014] The selected hard rock constitutive model comprises three parts: incremental calculation, strength failure criterion, and corresponding plastic stress correction, which are described in detail below:

[0015] (1) CWFS model

[0016] The CWFS model is based on the important characteristic of cohesion weakening and frictional strengthening during the brittle failure process of rock masses, and is proposed on the basis of the Mohr-Coulomb strength criterion. This model assumes that c and φ are functions of plastic strain, and the relationship between the values ​​of c and φ and plastic strain is shown in the following formula:

[0017] (1)

[0018] Where τ is the shear strength of the rock mass; σ n This represents the normal stress on the failure surface of the rock mass. j=1,3 represents the plastic strain increment; c and φ represent the cohesion and internal friction angle of the rock mass, respectively. This is the equivalent plastic strain.

[0019] (2)

[0020] in, , and The main plastic strain increment.

[0021] For ease of application in the FLAC3D numerical calculation software, referring to the FLAC3D user manual, the equivalent plastic strain calculation adopts the calculation formula recommended by the software:

[0022] (3)

[0023] (4)

[0024] (2) Brittle strength criterion

[0025] After the excavation of a high-stress tunnel, the tangential compressive stress σ1 at its sidewalls continuously increases. When σ1 exceeds a certain threshold, macroscopic cracks appear in the surrounding rock parallel to the direction of σ1. Figure 1 The failure mode of brittle fracture is called brittle failure, and the corresponding threshold is called "brittle strength." In engineering practice, the determination of brittle strength generally relies on empirical judgment or combined with laboratory testing methods, using prism-shaped rock samples with relatively small thicknesses to determine the brittle strength. Regression statistical analysis of a large amount of experimental data has revealed that, in the principal stress space, the influence of the middle principal stress on the rock strength conforms to the characteristics of a parabola. Therefore, its empirical brittle fracture strength satisfies the parabolic equation:

[0026] (5)

[0027] In the formula, σ1 is the brittle strength; σ2 is the intermediate principal stress; UCS is the average uniaxial compressive strength of the standard cylindrical specimen; K is the degradation factor of the rock in the field, with a value range of (0, 1], and 1 for indoor specimens. For the rock mass in the field, it is determined by back analysis of the measured data; parameters A, B, and C are determined by fitting the indoor biaxial brittle test data or by combining the suggested reference values ​​of existing studies.

[0028] Formula (6) yields the criterion for brittle strength:

[0029] (6)

[0030] When the stress state of a certain element is σ1 / σ3≥20, and When the value is greater than 0, brittle fracture will occur.

[0031] The above principles and methods are encapsulated according to the requirements of FLAC3D's constitutive model for use in FLAC3D numerical calculations.

[0032] It should be noted that establishing a new constitutive model in FLAC3D should follow the same computational rules as the original constitutive model, that is, calculating the stress increment at time step ∆t and the stress state update value at time t+∆t by using the stress state at time t and the total strain increment at time step ∆t. Therefore, the CWFS model and brittle constitutive model form of the tunnel hard and brittle rock mass under high ground stress are adopted.

[0033] The selected parameters of the hard rock brittle constitutive model, namely the rock mechanics parameters of the numerical model of the present invention, include: deformation parameters such as deformation modulus E and Poisson's ratio v, basic failure parameters such as cohesion c and internal friction angle φ, degradation factor K, parameters A, B, C, and brittle failure parameters such as uniaxial compressive strength UCS.

[0034] Step 2 detailed explanation:

[0035] The absolute error between the total number of computational grid cells within the measured brittle failure zone of the current excavation step of the deep-buried tunnel and the total number of computational grid cells entering the plastic state within the measured brittle failure zone after calculation by the FLAC3D numerical model is selected as the optimization objective function, which is explained in detail below:

[0036] By iteratively calculating and searching for the optimal parameter combination, the total number of computational grid elements entering the plastic state within the brittle fracture zone calculated by FLAC3D gradually approaches the total number of computational grid elements within the measured brittle fracture zone, thereby progressively correcting the parameters to be determined. The objective function can be set as follows:

[0037] (7)

[0038] In the formula: N1 and N2 are the total number of computational grid cells in the measured brittle fracture zone and the total number of computational grid cells that enter the plastic state in the measured brittle fracture zone after calculation by the FLAC3D numerical model, respectively.

[0039] Step 3 detailed explanation:

[0040] The selected SSA-IVM joint optimization algorithm is described in detail below:

[0041] (1) The SSA sparrow search optimization algorithm is a novel metaheuristic optimization algorithm inspired by the foraging and anti-predation behaviors of sparrows. Its principle is as follows:

[0042] Step 1: Initialize the population, number of iterations, and the ratio of predators to participants.

[0043] Step 2: Calculate fitness values ​​and sort them.

[0044] Step 3: Update the predator's location.

[0045] Step 4: Update the joiner's position.

[0046] Step 5: Update the location of the vigilant.

[0047] Step 6: Calculate fitness values ​​and update sparrow positions.

[0048] Step 7: Check if the stopping condition is met. If it is, exit and output the result. Otherwise, repeat Step 2-6.

[0049] The specific optimization rules for the SSA algorithm are as follows:

[0050] a) The discoverer typically possesses high energy reserves and is responsible for searching for areas with abundant food throughout the population, providing foraging areas and directions for all newcomers. In the model, the level of energy reserves depends on the fitness value of the individual sparrow.

[0051] b) Once a sparrow spots a predator, it begins to chirp as an alarm signal. When the alarm threshold exceeds the safe threshold, the sparrow will lead the newcomer to another safe area to forage.

[0052] c) The roles of discoverers and joiners are dynamic. Every sparrow can become a discoverer as long as it can find a better food source, but the proportion of discoverers and joiners in the entire population remains constant. In other words, if one sparrow becomes a discoverer, another sparrow will inevitably become a joiner.

[0053] d) The lower the energy level of the newcomers, the worse their foraging position within the population. Some hungry newcomers are more likely to fly to other places to forage for more energy.

[0054] e) During foraging, participants are always able to find the discoverer who provides the best food, and then obtain food from the best food or forage around that discoverer. At the same time, some participants may constantly monitor the discoverer in order to increase their predation rate and compete for food resources.

[0055] f) When they sense danger, sparrows on the edge of the group will quickly move to a safe area to get a better position, while sparrows in the middle of the group will move randomly to get closer to other sparrows.

[0056] In the simulation experiment, virtual sparrows are needed to find food. A population of n sparrows can be represented as follows:

[0057] (8)

[0058] Where d represents the dimension of the variable to be optimized, and n is the number of sparrows. Then, the fitness values ​​of all sparrows can be expressed in the following form:

[0059] (9)

[0060] Where f represents the fitness value.

[0061] In SSA (Specialized Social Association), the discoverer with a better fitness value will have priority in acquiring food during the search process. Furthermore, because the discoverer is responsible for finding food for the entire sparrow population and providing foraging directions for all participants, the discoverer has a larger foraging search area than the participants. According to rules (a) and (b), the discoverer's position update in each iteration is described as follows:

[0062] (10)

[0063] Where t represents the current iteration number, and j = 1, 2, 3, ..., d. max X is a constant representing the maximum number of iterations. ij This represents the position information of the i-th sparrow in the j-th dimension. α∈(0,1] is a random number.

[0064] R2 (R2∈[0,1]) and ST (ST∈[0.5,1]) represent the warning value and the safety value, respectively. Q is a random number that follows a normal distribution. L represents a 1×d matrix where every element is 1.

[0065] When R2 < ST, it means that there are no predators in the foraging environment, and the finder can perform extensive search operations.

[0066] If R2≥ST, it means that some sparrows in the population have discovered the predator and alerted the other sparrows in the population. At this time, all sparrows need to quickly fly to other safe places to forage.

[0067] For newcomers, rules (c) and (d) must be followed. As described earlier, during foraging, some newcomers will constantly monitor the discoverer. Once they realize that the discoverer has found better food, they will immediately leave their current position to compete for it. If they win, they can immediately obtain the discoverer's food; otherwise, they need to continue following rule (d). The location update of newcomers is described as follows:

[0068] (11)

[0069] Among them, X p This is currently the optimal position occupied by the discoverer, X. worst This represents the current worst position globally. A represents a 1×d matrix where each element is randomly assigned a value of 1 or -1, and A... + = A T (AA T ) -1 .

[0070] When i > n / 2, this indicates that the i-th entrant with a lower fitness value has not obtained food and is in a state of extreme hunger. At this time, it needs to fly to other places to find food in order to obtain more energy.

[0071] In the simulation, we assume that 10% to 20% of the sparrows are aware of the danger. Their initial locations are randomly generated within the population. According to rule (e), its mathematical expression can be represented as follows:

[0072] (12)

[0073] Among them, X best This is the current global optimal position. β, as the step size control parameter, is a random number following a normal distribution with a mean of 0 and a variance of 1. K∈[−1,1] is a random number, f i This is the fitness value of the current individual sparrow. g and f w These are the current best and worst fitness values ​​globally, respectively. ε is a constant to avoid zero in the denominator.

[0074] For simplicity, when f i >f g This indicates that the sparrows are currently on the fringes of their population and are extremely vulnerable to predators. X best This indicates that this location is the best and safest spot for the sparrows in the population. i =f g This indicates that sparrows in the middle of the population are aware of danger and need to move closer to other sparrows to minimize their risk of being preyed upon. K represents the direction of the sparrow's movement and is also a step size control parameter.

[0075] (2) The IVM algorithm is a fast sparse Gaussian process method based on information theory criteria. Its basic implementation steps are mainly divided into five steps: defining the marginal approximation function, calculating the posterior distribution, selecting the active subset, learning the kernel parameter θ, and data prediction. The specific principles are as follows:

[0076] a) Define the marginal approximation function

[0077] Based on Bayes' theorem, the joint distribution of the latent variable set f and the output observation y is expressed as:

[0078] (13)

[0079] In the formula, This is a noise model that illustrates the relationship between the latent variable f and the output observation y. Integrating the formula yields the marginal likelihood function:

[0080] (14)

[0081] In the formula, B is a diagonal matrix, and its nth diagonal element is , The value is .

[0082] A key assumption in Gaussian process models is that the noise model must follow a Gaussian distribution. Information Vector Machines (ADFs) construct approximate functions to replace non-Gaussian posterior distributions to maintain the model's applicability.

[0083] Therefore, the approximate formula for the marginal likelihood function of the regression problem is:

[0084] (15)

[0085] b) Calculate the posterior distribution

[0086] In the ADF algorithm, the process of approximating the true posterior distribution involves taking the training data n in J. i If added to I, the posterior distribution (f) can be updated as follows:

[0087] (16)

[0088] Minimizing the KL divergence using moment matching yields a new approximation: ,in yes The mean vector, yes The covariance matrices of these matrices are updated using the following formulas:

[0089] (17)

[0090] in, ,

[0091] Therefore, the ADF algorithm can be used to approximate any noise model using a Gaussian noise model. Least squares is a widely used model fitting method. For regression problems, this is equivalent to fitting a regression model with Gaussian noise, where each data point has a different inverse variance of the correlation.

[0092] (18)

[0093] Normalization constant Z n :

[0094] (19)

[0095] Log-normalization constant Approximate The parameter g in the formula for updating the mean vector and covariance matrix in and They can be represented as:

[0096] (20)

[0097] (twenty one)

[0098] Among them, g in G represents i The nth element, G in Then it means The nth diagonal element, then The nth diagonal element for

[0099] (twenty two)

[0100] In order to pass parameters and Approximate Gaussian distribution This can be summarized as likelihood, using parameters. Pair Gaussian target in By subdividing and rearranging the equations, we obtain:

[0101] (twenty three)

[0102] (twenty four)

[0103] c) Selection of activity subsets

[0104] In the information vector machine, the posterior differential entropy of each data point in J is calculated, and the data point with the maximum posterior differential entropy is selected and added to I. When the i-th information vector is selected for I, the posterior differential entropy of the n-th data point in J is:

[0105] (25)

[0106] For the formula To reduce memory usage, a sparse table is used, therefore the original variance matrix is ​​utilized. The continuous vector product can be solved to obtain :

[0107] (26)

[0108] In the formula, Let be an i×n matrix, with the k-th row being... , This represents the k-th information vector included in I. Because... covariance matrix The nth diagonal element in the matrix. Then the posterior covariance matrix. The update formula is expressed in diagonal form:

[0109] (27)

[0110] The formula for updating the posterior mean output vector is:

[0111] (28)

[0112] Therefore, the active subset can be determined using the information vector machine active subset selection algorithm:

[0113] a. Initialization. Set the number of information vectors in the activity subset to d; m=0; set μ=0; the inactive subset is all training sample data J; the active subset I is an empty set. It is an empty matrix.

[0114] b. When i=1, iterate through all training samples. Calculated according to formulas , , and . , and According to the formula Calculated .

[0115] c. Calculate the data points corresponding to the maximum a posteriori differential entropy. .

[0116] d. Based on formula and update and and using formula and calculation and .

[0117] e. Extension arrive Update again .

[0118] f. The nth i Add one data point to I and remove it from J.

[0119] g. Repeat steps a through f until d data points are selected and added to I.

[0120] By following the steps above, we can obtain the activity subset I used for training.

[0121] ④ Kernel parameter θ learning

[0122] The IVM algorithm uses training samples from the active subset to replace all training samples for Gaussian approximation, resulting in the marginal likelihood function.

[0123] (29)

[0124] The kernel function parameter θ is contained in In the middle, the scale conjugate gradient method is used to maximize the equation. The optimal θ can be obtained from the marginal likelihood function, as follows:

[0125] (30)

[0126] The training process of a Gaussian process involves directly selecting the kernel parameters from the training samples. Once the kernel parameters are determined, the samples can be predicted.

[0127] d) Data prediction

[0128] The optimal kernel parameters were calculated. Afterwards, according to The joint distribution of [ff(x)] can be obtained by integrating f with respect to the joint distribution and the posterior distribution, as follows:

[0129] (31)

[0130] in, It is the posterior covariance function;

[0131] This represents the posterior mean function. Once the distribution of f(x) is determined, the function value of f(x) at x can be predicted, which is the prediction process.

[0132] (3) The flowchart of the SSA-IVM joint optimization algorithm is shown below. Figure 3 The specific implementation steps are as follows:

[0133] Step A1: Set the SSA algorithm parameters: number of iterations T, population size N, predator P and entrant T; number of times to enter local optimization I. locad Maximum number of iterations I max The number of sparrows is NP.

[0134] Step A2: Set the IVM machine learning algorithm parameters: mean function d and covariance function (kernel function) in the active subset;

[0135] Step A3: Generate an initial random sparrow squad NP (i=1,2,3,…,NP), calculate the optimization objective function value f(X), sort the function values ​​in ascending order, and determine the current optimal position X of the food. L ;

[0136] Step A4: Perform global optimization using the SSA algorithm, X G X represents the optimal position obtained through a certain number of iterations. L This represents the optimal position obtained after the previous generation population update, i.e., the sparrow's position and its optimization objective function value, in the historical database X. record ;

[0137] Step A5: When the local iteration count of the SSA algorithm is I l Reaching the number of times to enter local optimization I locad Then, output X. L f(X) L ) and X record ;

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

[0139] Step A7: Train the IVM local proxy model using the training samples, that is, obtain the original optimization objective function value f(x) in X. best The distribution of the IVM approximate optimization objective function f in the local neighborhood IVM (x);

[0140] Step A8: Obtain the optimal value f in the distribution of the IVM approximate optimization objective function. IVM (x ib The corresponding x ib ;

[0141] Step A9: If f IVM (x ib It is better than f(X) L If X is updated, then X is updated. L= x ib ;

[0142] Step A10: Determine f(X) L If the convergence condition is met, the optimization ends; otherwise, return to step A4 until the global iteration count I = I max .

[0143] (4) The following steps are taken to identify the contour of the 3D brittle fracture zone using the SSA-IVM joint optimization algorithm combined with the feedback of the FLAC3D numerical model:

[0144] ① Select n model parameters from the numerical model parameters as optimization variables and determine their optimization range. Select other model parameters from their existing ranges. Determine the convergence condition of the optimization objective function f(X) as less than a threshold value f. Best ;

[0145] ② Based on the SSA-IVM joint optimization algorithm, determine the SSA algorithm parameters: number of iterations T, population size N, predator P and entrant T; number of times to enter a local optima I. locad Maximum number of iterations I max The number of sparrows is NP.

[0146] ③ Generate an initial random sparrow flock NP (i=1,2,3,…,NP), and perform FLAC3D forward calculation on the sparrow flock to obtain the optimization objective function value f(X). Sort the sparrows according to their fitness function values ​​from smallest to largest, and determine the current optimal position X of the food. L ;

[0147] ④ Use the SSA optimization algorithm for global optimization and record the sparrow's position information for each iteration, i.e., the sparrow's position and the objective function value calculated by FLAC3D, and store it in the historical database X. record ;

[0148] ⑤ When the local iteration count of the SSA algorithm is I l Reaching the number of times to enter local optimization I locad Then, output X. L f(X) L ) and X record ;

[0149] ⑥ Select X record Mid-range X L Location information of sparrows that are 3 to 5 times more numerous recently was used as training samples;

[0150] ⑦ Train the IVM local proxy model using the training samples, that is, obtain the original optimization objective function value f(x) in X. best The distribution of the IVM approximate optimization objective function f in the local neighborhood IVM (x);

[0151] ⑧ Obtain the objective function value through FLAC3D forward calculation, if f IVM (x ib It is better than f(X) L If X is updated, then X is updated. L= x ib ;

[0152] ⑨ If f(X) L ) < f BestIf the result is satisfactory, the optimization process ends and the optimal numerical model parameters are output; otherwise, return to step ④ until the iteration count I = I. max .

[0153] The selected SSA-IVM joint optimization algorithm, through I local optimization attempts... locad To terminate the SSA global optimization and enter the IVM local proxy model, I max with I locad The ratio should ideally be between 30 and 50.

[0154] Compared with the prior art, the beneficial effects of the present invention are:

[0155] (1) The IVM machine learning algorithm used in this invention to train the local surrogate model has the characteristics of adapting to small samples and having good regression and generalization. Due to the highly nonlinear mapping relationship between the mechanical parameters and mechanical response of rock mass, especially when the dimensionality of the mechanical parameters in the back analysis is high, using existing technology, i.e., using machine learning algorithms to train the global surrogate model, is difficult to guarantee a certain approximate accuracy under small sample conditions, which easily leads to optimization failure. Conversely, in order to improve the surrogate accuracy, a large number of training samples are generated, which also leads to the problem of excessive numerical reanalysis. Therefore, by using the IVM machine learning algorithm to train the local surrogate model, and performing surrogate in the local neighborhood of the objective function during the optimization process, the contradiction between numerical reanalysis and surrogate accuracy can be reduced.

[0156] (2) The SSA optimization algorithm used in this invention is used for global optimization of rock mechanics parameters inversion in deep-buried tunnels. It has the characteristics of high robustness, fast convergence speed and stronger global optimization ability. Compared with existing algorithms, the SSA optimization algorithm does not require complex parameter adjustment or operator design, is easy to implement and apply, and can find better surrounding rock mechanics parameters, making the difference between the measured brittle failure zone size parameters and the numerical simulation size parameters smaller, thus improving the reliability of rock mechanics parameters.

[0157] (3) The SSA-IVM joint optimization algorithm proposed in this invention, combined with the FLAC3D numerical model, is used for feedback identification of brittle fracture zone size parameters. The SSA optimization algorithm is used to globally optimize the brittle fracture zone size parameters, and FLAC3D calculates the fitness value. After a certain number of optimization iterations, sparrow information or historical sparrow information within a certain temperature range from the current optimal value in the historical database is selected to establish a training model and train a local surrogate model based on IVM. This obtains the IVM approximate fitness function distribution of the original fitness function in the local neighborhood, thereby finding a position better than the current optimal sparrow fitness. That is, the average difference between the brittle fracture zone size parameters measured in the current actual measurement and the brittle fracture zone size parameters obtained by FLAC3D numerical simulation is smaller, enabling it to quickly find the mechanical parameters of the brittle fracture zone rock mass that meet the convergence conditions. Compared with the existing single optimization algorithm or the use of a global surrogate model, this method effectively improves the efficiency of optimization back analysis while ensuring the calculation accuracy. Attached Figure Description

[0158] Figure 1 This is a schematic diagram showing the location of the brittle failure zone in a deep-buried tunnel engineering, provided by the present invention.

[0159] Figure 2 This is a schematic cross-sectional view of the brittle failure zone in a deep-buried tunnel engineering project provided by the present invention.

[0160] Figure 3 This is a flowchart of the SSA-IVM joint optimization algorithm provided in this invention.

[0161] Figure 4 This is a flowchart of an intelligent joint optimization and inversion method for mechanical parameters of surrounding rock of deep-buried tunnels based on three-dimensional brittle failure zone contours, provided by the present invention.

[0162] Figure 5 This is a schematic diagram of the simulation results of the brittle failure zone of the corresponding measuring point on the sidewall of a deep-buried tunnel project provided in Embodiment 2 of the present invention;

[0163] Figure 6 This is a two-dimensional model diagram of the IVM local optimization experience dataset provided in Embodiment 2 of the present invention;

[0164] Figure 7 This invention compares the number of FLAC3D calls and search time of two algorithms for searching the global optimal parameters of a deep-buried tunnel project. Detailed Implementation

[0165] The specific embodiments of the present invention will be further described and illustrated below with reference to the accompanying drawings and examples. It should be noted that the drawings only show the parts relevant to the present invention and not all results. Furthermore, the specific examples are only for explaining the present invention and not for limiting its scope.

[0166] Example 1

[0167] To understand the performance of the SSA-IVM joint optimization algorithm, Example 1 compares the proposed SSA-IVM joint optimization algorithm with the IVM optimization algorithm, the Genetic Algorithm (GA), the Particle Swarm Optimization (PSO) algorithm, and the Dung Beetle Optimization (DBO) algorithm based on single-peak and multi-peak benchmark functions. For all algorithms, the maximum number of iterations is 15000, and the number of specified operators is 50.

[0168] The benchmark function based on a single peak is:

[0169] (32)

[0170] The operator search range is [-200, 200], f min The value is 0, and f is used in this case. best Set it to 0.01. The single-peaked test function has only one optimum. After the algorithm completes the maximum number of iterations, the comparison algorithms all reach f. best The number of iterations is used to evaluate the convergence and exploratory capabilities of the algorithm.

[0171] The benchmark function based on multiple peaks is:

[0172] (34)

[0173] The operator search range is [-5.12, 5.12], and the smaller f2(x) is, the better. min The value is 0. The multi-peaked test function has multiple optimal values, but only one global optimum. The rest are local optima. After the algorithm completes the maximum number of iterations, the optimal values ​​reached by the algorithm are compared to evaluate the algorithm's ability to move from exploring local optima to reaching the global optimum.

[0174] Other parameters of the algorithms used in this comparison are shown in Table 1. For a detailed explanation of the optimization process of the SSA-IVM joint optimization algorithm, please refer to step 4. Figure 3 .

[0175] Table 2 shows the comparison of the number of calls to typical functions among the algorithms compared in this study. The SSA-IVM algorithm has a lower computational cost than the IVM, GA, PSO, and DBO algorithms, while exhibiting stronger global optimization capabilities. This indicates that the SSA-IVM joint optimization algorithm has significant advantages in typical test functions.

[0176] Table 1 Control parameter values ​​for each algorithm

[0177]

[0178] Table 2 Comparison of typical function call counts and optimal values ​​for each algorithm

[0179]

[0180] Example 2

[0181] Example 2 is a deep-buried tunnel project with a total length of 8.8km. The measured maximum principal stress in the tunnel sidewall area is about 30-32MPa, which is a high ground stress zone, and the brittle failure of the sidewall is obvious.

[0182] Figure 4 This is a schematic diagram of the simulation results of the brittle failure zone at the corresponding measuring point of the sidewall of a deep-buried tunnel project provided in an embodiment of the present invention. The specific method is as follows:

[0183] Step S1: Establish the FLAC3D numerical model of the tunnel. In this embodiment, based on geological survey data and tunnel design drawings, a computational mesh for the tunnel excavation body is established, such as... Figure 5 As shown. The computational grid for the tunnel region uses the hard rock CWFS constitutive model from step 1.

[0184] Step S2: Determine the optimization variables. Example 2: Based on the changes in the excavation location and the deformation and failure characteristics of the rock mass, rock mechanics parameters with significant variations are selected sequentially. Specifically, the degradation factor K, deformation modulus E, Poisson's ratio v, internal friction angle φ, and cohesion c are chosen as optimization variables. Among these, deformation modulus E, Poisson's ratio v, cohesion c, internal friction angle φ, and tensile strength σ are... t The optimal range was determined by referring to geological exploration data. The parameters of the brittle constitutive model of the tunnel sidewall, parameters A, B, and C, and the values ​​of the uniaxial compressive strength UCS were determined by on-site sampling and indoor testing, as shown in Table 3.

[0185] Table 3. Search range of parameters for the inversion model in the embodiment and other parameter values.

[0186]

[0187] Step S3: Establish the optimization objective function. The average difference between the brittle failure zone measured in the current excavation step and the brittle failure zone obtained from FLAC3D numerical simulation is used as the optimization objective function. Based on the specific instructions and formulas in Step 3, add an optimization objective function command flow to the FLAC3D numerical calculation script.

[0188] Step S4: Search for globally optimal parameters. When a tunnel excavation step is completed, based on the brittle failure zone of the sidewalls, the SSA-IVM joint optimization algorithm is used in conjunction with the tunnel FLAC3D numerical model to search for the optimal numerical model parameters.

[0189] Specifically, the global optimal parameter search based on the SSA-IVM-FLAC3D joint optimization algorithm was completed, with an optimized objective function value of 0.027. The results for the optimal numerical model parameters are shown in Table 4. A comparison of the number of FLAC3D calls and search time for the four algorithms is shown in [Table 4]. Figure 6 .

[0190] Table 4 Results of Parameters of Brittle Constitutive Model for Tunnel Sidewall in Examples

[0191]

[0192] Step S5: Based on the optimal model parameters of the tunnel FLAC3D numerical model, simulate the subsequent tunnel excavation process and obtain the predicted value of the brittle failure zone of the sidewall.

[0193] Specifically, after the global optimal parameter search is completed, the optimal model rock mechanics parameters are substituted into the tunnel FLAC3D numerical model to conduct numerical simulation of the next excavation cycle and obtain the predicted value of the brittle failure zone of the sidewall.

[0194] Example 2 describes a method for inverting the mechanical parameters of the surrounding rock in a tunnel. The specific implementation process is detailed in [link to example]. Figure 3 A novel intelligent algorithm based on Sparrow Optimization Algorithm (SSA) and Information Vector Machine (IVM) machine learning is proposed, combining the low-computational-cost IVM surrogate model with the powerful global optimization capability of SSA. Furthermore, this algorithm is integrated with the numerical computation software FLAC3D to propose an SSA-IVM-FLAC3D method for rock mass mechanics parameter inversion. The specific implementation steps are as follows:

[0195] Initialization algorithm parameter settings: SSA algorithm parameters: N=100, safety threshold ST=0.8, number of discoverers 50%, number of participants 30%, number of scouts 20%; After 7 iterations, the SSA algorithm starts IVM machine learning, selecting a total of 3×80=240 learning samples in the neighborhood of the locally optimal individual. Each time, IVM machine learning selects 100 samples from these 240 as information vectors, and repeats this selection process 8 times to obtain the 100 most "useful" information vectors for learning regression. The convergence criterion for the objective function of all three algorithms is ε=1×10⁻⁶. -3 ;

[0196] The range of the individual search space of the algorithm is set, which is the range of values ​​of the surrounding rock mechanical parameters to be inverted (see Table 5).

[0197] Table 5 Search range for surrounding rock mechanical parameters

[0198]

[0199] The optimization algorithm is started to perform optimization calculations. After an individual evolves according to the algorithm's evolution strategy, the individual's new position (80 sets of surrounding rock mechanical parameters) is stored in the data interface file A.

[0200] The FLAC 3D numerical calculation software is launched by calling a custom software command. The surrounding rock mechanical parameters in interface file A are read by the FISH program embedded in FLAC 3D and substituted into the established FLAC 3D numerical model to obtain the calculated brittle failure zone N2. N2 and the "measured three-dimensional dimensions of the brittle failure zone" N1 are then substituted into the objective function. The objective function values ​​corresponding to 80 sets of surrounding rock mechanical parameters were obtained and stored in data file B.

[0201] The objective function values ​​in data file B are read using a MATLAB program to obtain the fitness values ​​of all individuals. The minimum fitness value is selected by comparison, and the position coordinates of the individual corresponding to the minimum fitness value are the optimal combination of surrounding rock parameters for this generation of individuals.

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

[0203] The parameter inversion results of the two optimization algorithms are shown in Table 6. It can be seen that, compared with the SSA algorithm, the parameter inversion results of the SSA-IVM algorithm are closer to the true values.

[0204] Table 7 lists the cost and computation time of parameter inversion using different algorithms. As can be seen from Table 7, the computation cost and computation time of SSA-IVM are much lower than those of SSA, indicating that the SSA-IVM joint optimization algorithm has significant advantages in performing mechanical parameter inversion of surrounding rock of large cavern groups.

[0205] Table 6 Inversion results of surrounding rock mechanical parameters

[0206]

[0207] Table 7 Comparison of computation time between the two algorithms

[0208]

[0209] It should be noted that the purpose of disclosing the above examples is to help further understand the present invention. However, those skilled in the art will understand that various obvious changes, readjustments, and substitutions to the present invention will not depart from the scope of protection of the present invention. Therefore, the present invention is not limited to the content disclosed in the examples, and the scope of protection of the present invention is defined by the claims.

Claims

1. A method for inverting the mechanical parameters of surrounding rock in deep-buried tunnels based on a joint optimization algorithm of three-dimensional brittle failure zone profile and SSA-IVM, characterized in that, Includes the following steps: Step 1: Establish the FLAC3D numerical model of the tunnel.

2. Step 2: Establish the optimization objective function. The absolute error between the total number of computational grid cells within the measured brittle fracture zone and the total number of computational grid cells entering the plastic state within the measured brittle fracture zone after calculation by the FLAC3D numerical model is used as the optimization objective function.

3. Step 3: Search for the global optimum using an optimization algorithm. Using the tunnel surrounding rock mechanical parameters as optimization variables and aiming to achieve the global minimum value of the objective function, a joint optimization algorithm combining the Sparrow Optimization Algorithm (SSA) and Information Vector Machine (IVM), i.e., the SSA-IVM joint optimization algorithm, is used to perform global optimization, thereby obtaining the global optimum solution, i.e., the reasonable surrounding rock mechanical parameters.

4. The intelligent SSA-IVM joint optimization and inversion algorithm for mechanical parameters of surrounding rock in deep-buried tunnels based on brittle failure zone profiles as described in claim 1, characterized in that, The globally optimal search parameters are: Step A1: Set the SSA algorithm parameters: number of iterations T, population size N, predator P and entrant T; number of times to enter local optimization I. locad Maximum number of iterations I max The number of sparrows is NP. Step A2: Set the IVM machine learning algorithm parameters: mean function d and covariance function (kernel function) in the active subset; Step A3: Generate an initial random sparrow squad NP (i=1,2,3,…,NP), calculate the optimization objective function value f(X), sort the function values ​​in ascending order, and determine the current optimal position X of the food. L ; Step A4: Perform global optimization using the SSA algorithm, X G X represents the optimal position obtained through a certain number of iterations. L This represents the optimal position obtained after the previous generation population update, i.e., the sparrow's position and its optimization objective function value, in the historical database X. record ; Step A5: When the local iteration count of the SSA algorithm is I l Reaching the number of times to enter local optimization I locad Then, output X. L f(X) L ) and X record ; Step A6: Select X record Mid-range X L A certain range of sparrow information is used as training samples; Step A7: Train the IVM local proxy model using the training samples, that is, obtain the original optimization objective function value f(x) in X. best The distribution of the IVM approximate optimization objective function f in the local neighborhood IVM (x); Step A8: Obtain the optimal value f in the distribution of the IVM approximate optimization objective function. IVM (x ib The corresponding x ib ; Step A9: If f IVM (x ib It is better than f(X) L If X is updated, then X is updated. L= x ib ; Step A10: Determine f(X) L Check if the convergence condition is met; if so, end the optimization search. Otherwise, return to step A4 until the global iteration count I = I max .

5. The intelligent optimization and inversion algorithm for the mechanical parameters of surrounding rock in deep-buried tunnels based on the three-dimensional brittle failure zone profile as described in claim 1, characterized in that, The fitness function evaluation method is as follows: The optimal parameter combination is searched through iterative calculation. The absolute error between the total number of computational grid cells within the measured brittle fracture zone and the total number of computational grid cells entering the plastic state within the measured brittle fracture zone after calculation by the FLAC3D numerical model is used as the optimization objective function. The objective function can be set as follows:

6. In the formula: N1 and N2 are the total number of computational grid cells in the measured brittle fracture zone and the total number of computational grid cells that enter the plastic state in the measured brittle fracture zone after calculation by the FLAC3D numerical model, respectively.

7. The SSA-IVM-FLAC3D joint method according to claim 1, characterized in that, The training samples for the IVM local proxy model should preferably be selected from X. record Mid-range X L Location information for sparrows that are 3 to 5 times more numerous recently.

8. The SSA-IVM joint optimization algorithm according to claim 1, characterized in that, Through local optimization number I locad To stop the SSA global optimization and enter the IVM local proxy model, the selected I max with I locad The ratio should ideally be between 30 and 50.

Images