A method for TBM surrounding rock dynamic evolution modeling and parameter sensitivity analysis by combining numerical simulation and depth operator
Patent Information
- Application Number
- CN202611433884.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-09-16
- Publication Date
- 2026-10-23
AI Technical Summary
以有限元为代表的数值模拟方法能精细求解三维岩体应力-应变场，但单次非线性计算耗时长，难以支撑多参数组合、多开挖步工况下的批量分析需求
[0042]本发明实现了从高保真数据生成、物理约束代理建模、三维可视化评估、卡机风险评估、灾变因果机理提取到全局敏感性量化的完整技术方案，为TBM穿越断层破碎带时的围岩稳定性分析与掘进参数优化提供高效、高精度的科学分析手段。
Smart Images
Figure CN122965239A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of tunnel boring machine (TBM) construction technology, and in particular to a method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators. Background Technology
[0002] When tunnel boring machines (TBMs) traverse fault fracture zones and other adverse geological conditions, the stability of the surrounding rock becomes a significant issue. Accurately assessing the displacement field, stress field, and spatial distribution of the plastic zone of the surrounding rock under different excavation conditions and geological parameters, and identifying the key parameters affecting surrounding rock stability and their interaction effects, is crucial for ensuring the safe and efficient excavation of TBMs. Currently, related technologies and methods mainly suffer from the following shortcomings:
[0003] (1) The contradiction between the computational efficiency of numerical simulation and the expressive power of surrogate models. Numerical simulation methods, represented by the finite element method, can solve the stress-strain field of three-dimensional rock mass in detail, but the time-consuming single nonlinear calculation makes it difficult to support the batch analysis needs under multiple parameter combinations and multiple excavation steps. Existing surrogate model methods, such as response surface methodology and Kriging interpolation, can reduce computational costs, but their output is usually limited to a single scalar index, which cannot fully describe the spatial distribution characteristics of the three-dimensional physical field and cannot meet the visualization requirements of the global displacement field and stress field for surrounding rock stability analysis.
[0004] (2) Deep learning proxy models lack physical mechanism constraints. In recent years, deep learning methods such as deep operator networks have been introduced into the field of fast reasoning in physics. By learning the mapping relationship between input parameters and output physics, they can replace traditional numerical solvers. However, pure data-driven deep operator networks rely only on data fitting during training and lack the embedding of physical laws such as the governing equations of solid mechanics. This results in insufficient prediction accuracy in parameter space extrapolation or sparse data regions, and the learned mapping relationship lacks physical interpretability.
[0005] (3) Computational bottlenecks and dimensional limitations of parameter sensitivity analysis. Traditional control variable methods or local sensitivity analysis methods change individual parameters one by one and repeatedly run numerical simulations. The computational cost increases linearly with the number of parameters, and cannot quantify the nonlinear interaction effects between parameters. The Sobol global sensitivity analysis method can quantify the main effects and interaction effects of parameters through variance decomposition, but it requires a large number of samples to support the convergence of variance estimation. Directly relying on the finite element solver to generate samples leads to excessive computational cost, which limits its application in high-dimensional rock mass parameter spaces.
[0006] (4) The evolution mechanism of surrounding rock disaster lacks quantitative causal analysis methods. Existing methods for analyzing the mechanism of surrounding rock disaster are mostly limited to qualitative descriptions of displacement values and plastic zone ranges. They cannot quantitatively extract the causal evolution path of the surrounding rock state along the excavation steps during the TBM excavation process from the data. It is difficult to accurately reveal the direction, lag sequence and causal intensity of the causal relationship between various state variables. The understanding of the law of disaster evolution is still mainly based on empirical judgment.
[0007] (5) Insufficient spatial resolution for jamming risk assessment. Existing jamming risk assessment methods usually only compare the total frictional resistance with the rated thrust of the propulsion system as a whole, and give a single risk index. They cannot identify the differences in jamming risk at different spatial locations in the shield area, nor can they provide spatial distribution information of the surrounding rock safety factor at the same time, which restricts the spatial accuracy of jamming risk early warning and the pertinence of support design.
[0008] (6) Poor model adaptability under conditions of scarce field data. High-fidelity field monitoring data for adverse geological sections is difficult to obtain and costly. There is a systematic deviation between the surrogate model trained solely by numerical simulation and the actual field response. Existing methods lack an effective mechanism for calibrating numerical simulation-driven surrogate models using limited field monitoring data, resulting in insufficient applicability and prediction confidence of the model under complex real-world conditions. Summary of the Invention
[0009] The purpose of this invention is to provide a method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators, aiming to solve or improve at least one of the above-mentioned technical problems.
[0010] To achieve the above objectives, the present invention provides the following solution:
[0011] A method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators includes the following steps:
[0012] A three-dimensional steady-state finite element model is established, and Latin hypercube sampling and steady-state finite element solution are performed. The output is a high-fidelity design data table containing three-dimensional spatial coordinates, excavation face location, fault geometric parameters, rock mass mechanical parameters and corresponding displacement field and Mises stress field.
[0013] Construct a Fourier-enhanced DeepONet deep operator network and train it using the high-fidelity design data table;
[0014] The trained Fourier-enhanced DeepONet deep operator network is used to infer the arbitrary excavation length, and the global physical field is obtained. Then, the plastic zone is statistically analyzed and visualized in three dimensions.
[0015] Based on the global physical field obtained through reasoning, the shield contact pressure distribution is calculated by normal projection of stress tensor, the total frictional resistance is obtained by integration, and the jamming risk index is calculated to assess the jamming risk.
[0016] Using the trained Fourier-enhanced DeepONet deep operator network as the inference engine, Sobol global sensitivity analysis was performed using the high-fidelity design data table to quantify the main effects and interaction effects of each parameter.
[0017] Furthermore, the Fourier-enhanced DeepONet deep operator network consists of a branch network, a backbone network, a Fourier feature embedding module, and an intermediate physical computation module;
[0018] The branch network receives a one-dimensional parameter vector composed of rock mechanics parameters, fault geometry parameters and excavation face position, and outputs a branch coefficient vector through a fully connected layer.
[0019] The backbone network receives three-dimensional spatial coordinates, maps them to a high-dimensional random Fourier feature space through the Fourier feature embedding module, and then outputs a basis function vector through a fully connected layer. The final output is the inner product of the branch coefficient vector and the basis function vector, which contains three displacement field components and a Mises stress value.
[0020] The intermediate physics calculation module obtains the strain tensor by taking the partial derivatives of the three displacement field components with respect to the spatial coordinates. Combining the linear elastic constitutive relation, the stress tensor is derived from the strain tensor, and then the Mises equivalent stress based on the displacement gradient is calculated.
[0021] Furthermore, the Fourier-enhanced DeepONet deep operator network training employs a composite loss function, which is composed of a weighted sum of a data fidelity term, a spatial gradient regularization term, and a physical consistency constraint term.
[0022] Furthermore, the data fidelity item is:
[0023]
[0024] in, The number of samples in a single training batch. and To obtain the true values of the displacement field and Mises stress from the design data table, and The displacement field and Mises stress predictions are directly output by the network. These are the weighting coefficients for the stress term;
[0025] The spatial gradient regularization term is:
[0026]
[0027] in, The displacement gradient is obtained by automatically differentiating the predicted displacement field output by the network with respect to the input spatial coordinates. To determine the true values of the displacement gradient in the design data table;
[0028] The physical consistency constraints are:
[0029]
[0030] in, The Mises equivalent stress is derived from the displacement gradient. The Mises stress is directly output by the network;
[0031] The composite loss function is: ,in, and These are weighting coefficients, with values ranging from 0.1 to 0.5.
[0032] Furthermore, the trained Fourier-enhanced DeepONet deep operator network is used to infer the arbitrarily long excavation, and the global physical field is obtained. Then, plastic zone statistics and 3D visualization are performed specifically as follows:
[0033] Based on the global stress field obtained through reasoning, the plastic state of each spatial sampling point is determined point by point using the Mohr-Coulomb yield criterion. The total volume of the plastic zone is obtained by summing the volumes of the Voronoi elements corresponding to all plastic sampling points. The minimum distance from the excavation face boundary to the farthest boundary of the plastic zone is determined by radial search as the loosening zone thickness. The displacement field, stress field, and plastic zone identifier are rendered by color mapping, and the distribution of surrounding rock displacement, stress concentration area, and spatial morphology of the plastic zone are displayed synchronously in the same three-dimensional view.
[0034] Furthermore, the risk assessment for the card reader specifically includes:
[0035] The stress tensor of the surrounding rock is extracted at each sampling point on the outer surface of the shield, and the contact pressure is obtained by calculating the projection of the stress tensor onto the outer normal direction of the shield surface. The contact pressure at all sampling points on the outer surface of the shield is integrated along the shield surface, and the total frictional resistance is obtained by combining the friction coefficient between the shield and the surrounding rock. The ratio of the total frictional resistance to the rated thrust of the TBM propulsion system is used as the jamming risk index. When the index exceeds the preset threshold, it is determined that there is a jamming risk at the current excavation length.
[0036] Furthermore, the Sobol global sensitivity analysis specifically includes:
[0037] Batch reasoning was performed using all parameter combinations in the high-fidelity design data table to extract the loosening zone thickness and plastic zone volume as catastrophic response indicators. The Sobol variance decomposition method was used to calculate the first-order sensitivity index and the total-order sensitivity index of each input parameter, quantifying the main effect of the individual parameter and the interaction effect between parameters.
[0038] Furthermore, it also includes the step of extracting the spatial causal evolution path of the excavation sequence of the surrounding rock disaster. Specifically, the global displacement field and stress field under each excavation length are obtained by reasoning through the trained Fourier-enhanced DeepONet deep operator network. The corresponding plastic zone range and loose zone thickness are statistically analyzed by combining the Mohr-Coulomb yield criterion, and the surrounding rock state evolution sequence is constructed according to the excavation sequence.
[0039] A 3D convolutional visual Transformer encoder is used to encode the spatial features of the three-dimensional physical field at each excavation length, obtain the corresponding high-dimensional feature vectors, and arrange them according to the excavation sequence. The excavation step sequence attention matrix is calculated to identify the key turning excavation lengths in the evolution of the surrounding rock state.
[0040] The PCMCI causal discovery algorithm is run with the excavation step sequence attention matrix as a prior constraint. False causal edges caused by common dependent variables are eliminated through conditional independence tests. A directed acyclic causal graph containing key state variables and directed causal edges is output. The causal evolution path of the surrounding rock disaster along the excavation step sequence is quantitatively extracted from the directed acyclic causal graph.
[0041] The present invention discloses the following technical effects:
[0042] This invention realizes a complete technical solution from high-fidelity data generation, physical constraint proxy modeling, three-dimensional visualization evaluation, machine jamming risk assessment, disaster causal mechanism extraction to global sensitivity quantification, providing efficient and high-precision scientific analysis methods for the surrounding rock stability analysis and tunneling parameter optimization when TBMs cross fault fracture zones. Attached Figure Description
[0043] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0044] Figure 1 This is an overall flowchart of the method proposed in the embodiments of the present invention;
[0045] Figure 2This is a diagram of the Fourier-DeepONet network architecture in an embodiment of the present invention;
[0046] Figure 3 This is a schematic diagram of the composite loss function structure in an embodiment of the present invention;
[0047] Figure 4 This is an example diagram of a 3D TBM tunnel visualization cloud map in an embodiment of the present invention;
[0048] Figure 5 This is a schematic diagram of the Sobol sensitivity index in an embodiment of the present invention. Detailed Implementation
[0049] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0050] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.
[0051] like Figures 1-5 As shown, this invention provides a method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators, including:
[0052] 1. Generate high-fidelity design data tables.
[0053] This step utilizes COMSOL Multiphysics to establish a three-dimensional steady-state finite element model of the fault fracture zone, and trains the research nodes through the built-in surrogate model to perform Latin hypercube sampling and batch solving, generating a high-fidelity design data table covering multiple parameter combinations and multiple excavation lengths.
[0054] (1) Establishment of a three-dimensional refined numerical model.
[0055] ① Geometry and Materials.
[0056] A three-dimensional geometric model was created in COMSOL, including the tunnel excavation face, surrounding rock, and fault fracture zone. The model's extent was set to five times the tunnel diameter outward from the tunnel center to eliminate boundary effects. The fault fracture zone was a plate-like region with a specific dip angle, strike, and width.
[0057] ② Physical interface and constitutive model.
[0058] The COMSOL solid mechanics physics interface was selected as the solution framework. For the fault fracture zone, a Rocks sub-node was added under the linear elastic material node, and the generalized Hoek-Brown yield criterion was used to describe its mechanical behavior. This model reflects the strength characteristics of the jointed rock mass and the structural plane-controlled failure mode through geological strength indices, intact rock parameters, and perturbation factors. For the relatively intact surrounding rock, a Soil Plasticity sub-node was added under the linear elastic material node, and the Mohr-Coulomb yield criterion was used for simulation.
[0059] ③ Grid division.
[0060] The computational domain uniformly adopts an unstructured tetrahedral mesh. Adaptive mesh refinement is implemented near the tunnel excavation face and within the fault fracture zone, with the minimum mesh size controlled at 0.1m and the mesh size relaxed to 1.0m in areas far away, to ensure accurate capture of stress concentration areas and strain localization phenomena.
[0061] ④ Boundary conditions and steady-state studies.
[0062] The model is fully constrained at the bottom, normally constrained on all four sides, and the self-weight stress of the overburden is applied at the top. The initial in-situ stress field is based on the lateral pressure coefficient applied according to measured in-situ stress data. A steady-state approach is adopted, with each excavation step corresponding to a specific tunnel excavation length, and an independent steady-state solution is performed at that length. The excavation step length is set according to the actual project progress.
[0063] (2) Proxy model setup and training.
[0064] ① Add a node for training and researching proxy models.
[0065] Create a new proxy model training study in the COMSOL model developer. This study node automatically creates and manages DOE (Design of Experiments) jobs at the underlying level.
[0066] ② Define the input parameters and their value range.
[0067] In the surrogate model training research settings window, specify the input parameters for sampling and their value ranges. Input parameters include: fault dip angle (60°–80°), fracture zone length (30–50 m), elastic modulus E (1–30 GPa), cohesion c (0.1–2.0 MPa), internal friction angle φ (15°–45°), and permeability coefficient k (…). ), geostress lateral pressure coefficient (0.5~2.5).
[0068] ③ Define the output focus.
[0069] Add output parameters of interest: total displacement (solid.disp) and Mises equivalent stress (solid.misesGp). These physical quantities will be recorded at all spatial sampling points for each parameter combination and each excavation length.
[0070] ④ Perform calculations and export the design data table.
[0071] Click "Compute," and COMSOL will automatically update the model parameters according to the LHS sampling scheme and perform steady-state finite element analysis for each excavation length. After the solution is complete, all input parameter combinations, spatial coordinates (x, y, z), and corresponding displacements and Mises stress values will be automatically summarized in the Design Data table. Export this table as a CSV file.
[0072] (3) Data preprocessing.
[0073] The exported design data table is cleaned, and rows with null values are removed. The input parameters (including coordinates) and output physical quantities are normalized to the interval [-1,1] or [0,1], respectively. The data is then randomly divided into training, validation, and test sets in a ratio of 7:2:1.
[0074] 2. Construction and training of the Fourier-DeepONet proxy model.
[0075] Using the high-fidelity design data table generated in step 1, a Fourier-enhanced DeepONet deep operator network with embedded physical constraints is trained to achieve millisecond-level inference of the displacement and stress fields of the surrounding rock.
[0076] (1) Network architecture design.
[0077] ① Branch Net.
[0078] Input: Rock mass mechanics parameters (E, c, φ, k, The fault geometry parameters (dip angle, fracture zone length) and the current position of the excavation face are concatenated into a one-dimensional vector. Network structure: 6-layer fully connected network, 512 neurons per layer, using the Swish activation function. Output: P-dimensional branch coefficient vector β(P), where P is the input parameter vector.
[0079] ② Trunk Net and Fourier feature embedding.
[0080] The backbone network receives three-dimensional spatial coordinates Introduce random Fourier feature maps in its first layer:
[0081]
[0082] in, Let B be the input three-dimensional spatial coordinate vector, where B ∈ R. m×3 The elements are from a Gaussian distribution Random sampling is used, with σ set to 1.0; m is the dimension of the Fourier feature map. The mapped high-dimensional vector is fed into a subsequent 5-layer fully connected network, with 512 neurons per layer. The activation function is GELU, and the output is a P-dimensional basis function vector τ. .
[0083] ③ Output fusion.
[0084] The final output of the network is the inner product of the branch coefficients and the backbone basis functions:
[0085]
[0086] in Includes three displacement components A Mises stress P is 256.
[0087] (2) Intermediate physical calculation module.
[0088] An explicit physics calculation module is embedded during the network forward propagation process, and the displacement-stress consistency relationship is derived using automatic differentiation.
[0089] ① Calculate strain using displacement gradient.
[0090] Using PyTorch's automatic differentiation function, spatial partial derivatives of the displacement components are calculated to obtain the strain tensor components:
[0091]
[0092] in, These are the i-th and j-th components of the small strain tensor calculated from the displacement field predicted by the network; and These are the components of the predicted displacement vector in the i-th and j-th coordinate directions, respectively; and These are spatial coordinate components; Represents displacement components Spatial coordinates The partial derivatives; i, j = 1, 2, 3, corresponding to the x, y, and z directions in the three-dimensional Cartesian coordinate system, respectively. When i = j, Indicates normal strain; when i ≠ j, This represents the tensor shear strain, and its corresponding engineering shear strain is... .
[0093] ② Determine stress using constitutive relations.
[0094] Calculate the stress tensor based on the linear elastic constitutive relation:
[0095]
[0096] in, For volumetric strain, The symbol for Kronecker. λ and μ are Lamé constants:
[0097]
[0098] Where ν is Poisson's ratio, taken as 0.25. This is the equivalent elastic modulus of the rock mass (Young's modulus).
[0099] ③ Derivation of Mises equivalent stress.
[0100] Calculation of Mises equivalent stress based on displacement gradient from stress tensor:
[0101]
[0102] (3) Design of composite loss function.
[0103] ① Data fidelity item.
[0104]
[0105] Where α = 1.0.
[0106] ② Spatial gradient regularization term.
[0107]
[0108] The gradient is obtained through automatic differentiation.
[0109] ③ Physical consistency constraints.
[0110]
[0111] ④ Weighted total loss.
[0112] The composite loss function is:
[0113]
[0114] in, =0.3, =0.3.
[0115] (4) Training configuration and convergence criteria.
[0116] Using the AdamW optimizer, initial learning rate A cosine annealing decay strategy is employed. Each batch contains 128 samples, and the total number of training rounds is [number missing]. When the validation set is relatively Training is stopped and the model weights are saved when the error does not decrease for 5000 consecutive rounds and falls below 2%.
[0117] 3. Post-processing statistics of the plastic zone and generation of 3D visualization cloud maps.
[0118] The trained surrogate model is used to infer the excavation length, and the global physical field is obtained. Then, the plastic zone is statistically analyzed and visualized in three dimensions.
[0119] (1) Mohr-Coulomb plasticity determination.
[0120] For each spatial sampling point, the maximum principal stress is calculated from the stress tensor output by the surrogate model. and minimum principal stress Substitute the yield function:
[0121]
[0122] when When the value is greater than 0, the point is marked as being in a plastic state.
[0123] (2) Calculation of the volume of the plastic zone and the thickness of the loose zone.
[0124] By summing the Voronoi units represented by each plastic point, the total volume of the plastic region is obtained. Within the tunnel cross-section, the plastic point is searched radially outward from the excavation face; the furthest distance is the thickness of the loosened zone. .
[0125] (3) Generation of 3D visualization cloud map.
[0126] The displacement amplitude is rendered as a blue-green-red gradient cloud map, the Mises stress is rendered as a blue-yellow-red gradient cloud map, and the plastic zone is overlaid with a semi-transparent red layer. The tunnel excavation face outline is displayed in the same 3D view, forming a composite cloud map of displacement, stress, and plastic zone.
[0127] 4. Card risk assessment and security factor distribution cloud map generation.
[0128] (1) Calculation of contact pressure.
[0129] Based on the surrounding rock stress tensor field output by the proxy model, stress tensors are extracted at various spatial sampling points on the outer surface of the shield. Calculate the stress tensor at that point in the direction of the outer normal to the shield surface. The contact pressure is obtained by projecting the image onto the surface.
[0130]
[0131] in, It is the unit outward normal vector, determined by the shield geometry.
[0132] (2) Calculation of total friction resistance and jamming risk index.
[0133] Integrate the contact pressure at all sampling points on the outer surface of the shield along the shield surface to obtain the total frictional resistance:
[0134]
[0135] in, The friction coefficient between the shield and the surrounding rock is taken as 0.3, and M is the number of sampling points on the shield surface. Let be the outer surface area of the shield represented by the i-th sampling point.
[0136] Card jam risk index:
[0137]
[0138] in, For the rated thrust of the TBM propulsion system, when When the threshold is exceeded, a risk of system crash is identified, and a safety factor distribution cloud map is generated.
[0139] (3) Overlay of comprehensive risk cloud map.
[0140] The distribution of the plastic zone, the cloud map of the location of the jamming risk, and the cloud map of the safety factor distribution are all superimposed on the three-dimensional tunnel model to form a comprehensive risk cloud map.
[0141] 5. Extraction of spatial causal evolution path of excavation steps in surrounding rock disaster.
[0142] (1) Generation of high-density excavation step sequence state sequence.
[0143] For critical sections where the TBM tunnels traverse fault fracture zones, such as the area from 5m from the front boundary of the fault to 5m from the rear boundary after exiting the fault, the current position of the excavation face is gradually changed in increments of 0.2 to 0.5m. For each excavation length, the corresponding rock mechanics parameters (elastic modulus E, cohesion c, internal friction angle φ, permeability coefficient k, and in-situ stress lateral pressure coefficient) are updated. The fault geometry parameters (dip angle, fracture zone length) and the current excavation face position are concatenated into an input vector, which is then input into the surrogate model trained and converged in step 2 for forward inference. This yields the displacement field (three components) and Mises stress field of all spatial sampling points within the entire three-dimensional computational domain at that excavation length. Simultaneously, using the Mohr-Coulomb yield criterion and plastic zone post-processing method described in step 3, the range, volume, and loosened zone thickness of the plastic zone at each excavation length are statistically analyzed and recorded as additional state information along with the displacement and stress fields. The aforementioned state information for all excavation lengths is arranged in the excavation sequence to form a sequence of surrounding rock state evolution.
[0144] (2) 3D convolutional visual Transformer encoding.
[0145] For the three-dimensional physical field corresponding to each excavation length in the above state evolution sequence, a 3D convolutional visual Transformer encoder is used for spatial feature compression. The encoder consists of two stacked parts: the first part is a three-layer 3D convolutional neural network, each layer using a 3×3×3 three-dimensional convolutional kernel with a stride of 1. Zero padding is applied to the input boundaries to maintain the spatial dimensions, and ReLU is used as the activation function. The number of channels in the three convolutional layers is set to 32, 64, and 128 respectively, extracting detailed features such as local stress concentration areas, displacement gradient abrupt change zones, and plastic zone boundaries in the three-dimensional displacement and stress fields layer by layer. The second part is a Transformer self-attention layer containing eight attention heads, each with a query, key, and value vector dimension of 64. The feature map output from the third 3D convolutional layer is flattened into a sequence and input into the self-attention layer. The multi-head self-attention mechanism captures the global spatial distribution features of the loosening zone and plastic zone in the entire three-dimensional computational domain, as well as the long-distance dependencies between regions. The encoder finally outputs a high-dimensional feature vector with a fixed dimension of 512 for the three-dimensional physical field at each excavation length. Arrange the feature vectors corresponding to all excavation lengths in the order of excavation to form a feature sequence of the evolution of the surrounding rock state to be analyzed.
[0146] (3) Calculation of the attention matrix for the excavation sequence.
[0147] For the sequence of surrounding rock state evolution characteristics obtained in step (2), calculate its self-attention matrix. Then, calculate the feature vector corresponding to the t-th excavation length. As the query vector, the feature vector corresponding to the τth excavation length As a key vector, the attention weights between the two are calculated by scaling the dot product:
[0148]
[0149] Where d=512 is the dimension of the feature vector, and T is the total excavation length. Attention weights The dependence of the surrounding rock state at the t-th excavation length on the surrounding rock state at the τ-th excavation length was quantified; a larger value indicates a stronger correlation between the states at the two excavation lengths. By traversing all t and τ, the self-attention matrix A for the entire T×T is obtained. This matrix can reveal key turning points in the evolution of the surrounding rock state—when the attention weight undergoes a sudden change or forms a local peak at a certain excavation length, it signifies a significant change in the surrounding rock state at that point, potentially a critical node in the catastrophic evolution.
[0150] (4) PCMCI causal discovery.
[0151] Using a self-attention matrix as the prior constraint graph, the PCMCI causal discovery algorithm is run. False causal edges caused by codependent variables are eliminated through a conditional independence test, with a significance level set at 0.05. Attention weights are used to prune candidate causal edges, retaining only causal associations with weights greater than a preset threshold. The final output is a directed acyclic causal graph containing nodes and directed edges. Nodes are automatically identified by the algorithm from the state evolution sequence, representing key state features in the evolution of the surrounding rock disaster. Directed edges are labeled with the number of lagging excavation steps and the causal strength coefficient. The complete evolution path of the surrounding rock disaster along the excavation steps can be quantitatively extracted from the causal graph.
[0152] 6. Global sensitivity analysis of Sobol.
[0153] (1) Use the existing design data table to perform batch inference of the surrogate model.
[0154] Input all parameter combinations (without adding new sampling points) from the design data table generated in step 1 into the trained surrogate model, and perform batch inference to obtain the displacement field and stress field corresponding to each sample point.
[0155] (2) Extraction of disaster response indicators.
[0156] Post-processing based on the yield criterion was performed to statistically analyze the loosened zone thickness for each parameter combination and each excavation length. and plastic zone volume , as an output variable of the sensitivity analysis.
[0157] (3) Sobol variance decomposition calculation.
[0158] Calculate each input parameter based on the samples contained in the design data table. First-order sensitivity Total order sensitivity :
[0159]
[0160] in, Disaster response indicators Total variance Quantization parameters The proportion of each individual effect to the total variance. Includes parameters With all other parameters The total contribution of the interaction effect; by comparing the parameters. and The numerical values are used to obtain a global ranking of the degree of influence of each parameter on the loosening zone thickness and the plastic zone volume, and to identify the main control parameters and parameter pairs with significant interaction effects. All expectation and variance calculations are based on all parameter combination samples naturally contained in the design data table and their corresponding surrogate model inference results.
[0161] (4) Sensitivity index analysis and master control parameter identification.
[0162] according to and The parameters with the greatest impact on the loose zone and plastic zone were determined by sorting their sizes, and the interaction effects between the parameters were identified. The analysis results are presented in tabular and bar chart formats.
[0163] The innovation and advantages of this invention:
[0164] 1. Fourier-DeepONet proxy model with embedded physical constraints.
[0165] A stochastic Fourier feature map is introduced into the input layer of the backbone network to enhance the ability to capture high-frequency localized features of displacement and stress fields. An explicit physical calculation module is embedded in the middle layer of the network. The strain tensor and stress tensor are derived from the displacement field by automatic differentiation, and then the Mises equivalent stress is calculated. The deviation between the derived stress and the stress directly output by the network is added to the loss function as a physical consistency constraint. Together with the data fidelity term and the spatial gradient regularization term, they form a composite loss function, which enables the surrogate model to have higher prediction accuracy and physical interpretability in sparse data regions and when extrapolating parameters.
[0166] 2. A surrounding rock stability assessment system that includes statistics on the plastic zone and three-dimensional visualization.
[0167] By combining the Mohr-Coulomb yield criterion to determine the plastic state point by point, and statistically analyzing the volume of the plastic zone and the thickness of the loosened zone, the distribution of the plastic zone is superimposed with the displacement cloud map and stress cloud map to generate a three-dimensional visualization cloud map, thereby realizing the synchronous three-dimensional display of surrounding rock deformation, stress concentration and failure area.
[0168] 3. Catastrophic evolution path extraction based on attention-enhanced causal discovery.
[0169] A high-density excavation step state sequence is generated using a proxy model as a high-speed inference engine. A 3D convolutional visual Transformer encoder is used for spatial feature compression. The excavation step attention matrix is calculated and used as a constraint to run the PCMCI causal discovery algorithm. The key state nodes of the surrounding rock disaster are automatically identified from the data, and the causal relationship direction and lag step number are quantitatively extracted.
[0170] 4. Global sensitivity analysis of Sobol based on the design data table.
[0171] Sobol variance decomposition was performed directly using the existing design data table generated by COMSOL, calculating the first-order and total-order sensitivity indices of each parameter, quantifying main effects and interaction effects, and identifying the main control parameters. The entire process did not require additional calls to the finite element solver, and the sample size of the design data table ensured the numerical stability of the variance decomposition.
[0172] This invention addresses the problems of low computational efficiency of traditional finite element methods, lack of physical constraints in proxy models, reliance on engineering simplification assumptions for jam risk assessment, and high computational cost of parameter sensitivity analysis when TBMs traverse fault fracture zones. It proposes a method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators. This method trains research nodes using a surrogate model built into COMSOL to perform Latin hypercube sampling, generating a high-fidelity design data table covering multiple parameter combinations and excavation lengths. Based on this, a deep operator network with embedded stochastic Fourier feature maps and intermediate physical calculation modules is constructed. Data fidelity, spatial gradient regularization, and physical consistency constraints are simultaneously applied to the loss function, training a surrogate model capable of inferring the three-dimensional displacement and stress fields of the surrounding rock at the millisecond level. The Mohr-Coulomb yield criterion is combined to perform post-processing statistics of the plastic zone and generate a three-dimensional visualization cloud map. The shield contact pressure is directly extracted based on the stress tensor normal projection to calculate the jamming risk index. A 3D convolutional visual Transformer encoder and attention-enhanced PCMCI causal discovery algorithm are used to quantitatively extract the surrounding rock disaster evolution path from the high-density excavation step sequence. Finally, existing design data tables are directly used for Sobol global sensitivity analysis to quantify the main effects and interaction effects of each parameter.
[0173] Compared to the traditional numerical simulation approach that repeatedly calls the finite element solver, this invention achieves efficient reasoning and comprehensive visualization assessment of the three-dimensional physical field of the surrounding rock through a surrogate model with embedded physical constraints. It directly quantifies the risk of tunnel boring machine jamming through a solid mechanics stress projection method, elevates the catastrophic mechanism of the surrounding rock from a qualitative description to a quantitative causal path expression along the excavation sequence through a causal discovery algorithm, and avoids additional finite element calculation costs by utilizing existing design data tables for sensitivity analysis. The resulting complete closed-loop technical solution, encompassing high-fidelity data generation, physical constraint surrogate modeling, three-dimensional visualization of the plastic zone, jamming risk assessment, catastrophic causal path extraction, and global sensitivity analysis, provides an efficient, high-precision, and highly interpretable full-process scientific analysis method for the stability analysis of surrounding rock and the optimization of tunneling parameters when TBMs traverse fault fracture zones.
[0174] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.
[0175] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.
Claims
1. A method for dynamic evolution modeling and parameter sensitivity analysis of TBM surrounding rock that integrates numerical simulation and depth operators, characterized in that, Includes the following steps: A three-dimensional steady-state finite element model is established, and Latin hypercube sampling and steady-state finite element solution are performed. The output is a high-fidelity design data table containing three-dimensional spatial coordinates, excavation face location, fault geometric parameters, rock mass mechanical parameters and corresponding displacement field and Mises stress field. Construct a Fourier-enhanced DeepONet deep operator network and train it using the high-fidelity design data table; The trained Fourier-enhanced DeepONet deep operator network is used to infer the excavation length, and the global physical field is obtained. Then, the plastic zone is statistically analyzed and visualized in three dimensions. Based on the global physical field obtained through reasoning, the shield contact pressure distribution is calculated by normal projection of stress tensor, the total frictional resistance is obtained by integration, and the jamming risk index is calculated to assess the jamming risk. Using the trained Fourier-enhanced DeepONet deep operator network as the inference engine, Sobol global sensitivity analysis was performed using the high-fidelity design data table to quantify the main effects and interaction effects of each parameter.
2. The TBM-based dynamic evolution modeling and parameter sensitivity analysis method for surrounding rock, which integrates numerical simulation and depth operators, as described in claim 1, is characterized in that... The Fourier-enhanced DeepONet deep operator network consists of a branch network, a backbone network, a Fourier feature embedding module, and an intermediate physical computation module; The branch network receives a one-dimensional parameter vector composed of rock mechanics parameters, fault geometry parameters and excavation face position, and outputs a branch coefficient vector through a fully connected layer. The backbone network receives three-dimensional spatial coordinates, maps them to a high-dimensional random Fourier feature space through the Fourier feature embedding module, and then outputs a basis function vector through a fully connected layer. The final output is the inner product of the branch coefficient vector and the basis function vector, which contains three displacement field components and a Mises stress value. The intermediate physics calculation module obtains the strain tensor by taking the partial derivatives of the three displacement field components with respect to the spatial coordinates. Combining the linear elastic constitutive relation, the stress tensor is derived from the strain tensor, and then the Mises equivalent stress based on the displacement gradient is calculated.
3. The TBM-based dynamic evolution modeling and parameter sensitivity analysis method for surrounding rock, which integrates numerical simulation and depth operators, as described in claim 1, is characterized in that... The Fourier-enhanced DeepONet deep operator network is trained using a composite loss function, which consists of a weighted sum of a data fidelity term, a spatial gradient regularization term, and a physical consistency constraint term.
4. The TBM-based modeling and parameter sensitivity analysis method for surrounding rock dynamics that integrates numerical simulation and depth operators as described in claim 3, is characterized in that... The data fidelity item is: in, The number of samples in a single training batch. and To obtain the true values of the displacement field and Mises stress from the design data table, and The displacement field and Mises stress predictions are directly output by the network. These are the weighting coefficients for the stress term; The spatial gradient regularization term is: in, The displacement gradient is obtained by automatically differentiating the predicted displacement field output by the network with respect to the input spatial coordinates. To determine the true values of the displacement gradient in the design data table; The physical consistency constraints are: in, The Mises equivalent stress is derived from the displacement gradient. The Mises stress is directly output by the network; The composite loss function is: in, and These are weighting coefficients, with values ranging from 0.1 to 0.
5.
5. The TBM-based dynamic evolution modeling and parameter sensitivity analysis method for surrounding rock that integrates numerical simulation and depth operators according to claim 1, characterized in that, The trained Fourier-enhanced DeepONet deep operator network is used to infer the global physical field for arbitrary excavation lengths. Then, plastic zone statistics and 3D visualization are performed as follows: Based on the global stress field obtained through reasoning, the plastic state of each spatial sampling point is determined point by point using the Mohr-Coulomb yield criterion. The total volume of the plastic zone is obtained by summing the volumes of the Voronoi elements corresponding to all plastic sampling points. The minimum distance from the excavation face boundary to the farthest boundary of the plastic zone is determined by radial search as the loosening zone thickness. The displacement field, stress field, and plastic zone identifier are rendered by color mapping, and the distribution of surrounding rock displacement, stress concentration area, and spatial morphology of the plastic zone are displayed synchronously in the same three-dimensional view.
6. The TBM-based dynamic evolution modeling and parameter sensitivity analysis method for surrounding rock, which integrates numerical simulation and depth operators, as described in claim 1, is characterized in that... The specific risks of the card machine assessment are as follows: The stress tensor of the surrounding rock is extracted at each sampling point on the outer surface of the shield, and the contact pressure is obtained by calculating the projection of the stress tensor onto the outer normal direction of the shield surface. The contact pressure at all sampling points on the outer surface of the shield is integrated along the shield surface, and the total frictional resistance is obtained by combining the friction coefficient between the shield and the surrounding rock. The ratio of the total frictional resistance to the rated thrust of the TBM propulsion system is used as the jamming risk index. When the index exceeds the preset threshold, it is determined that there is a jamming risk at the current excavation length.
7. The TBM-based dynamic evolution modeling and parameter sensitivity analysis method for surrounding rock, which integrates numerical simulation and depth operators according to claim 1, is characterized in that... The Sobol global sensitivity analysis specifically includes: Batch reasoning was performed using all parameter combinations in the high-fidelity design data table to extract the loosening zone thickness and plastic zone volume as catastrophic response indicators. The Sobol variance decomposition method was used to calculate the first-order sensitivity index and the total-order sensitivity index of each input parameter, quantifying the main effect of the individual parameter and the interaction effect between parameters.
8. The TBM surrounding rock dynamic evolution modeling and parameter sensitivity analysis method integrating numerical simulation and depth operator as described in claim 1, characterized in that, It also includes the step of extracting the spatial causal evolution path of the excavation sequence of the surrounding rock disaster. Specifically, the global displacement field and stress field under each excavation length are obtained by reasoning through the trained Fourier-enhanced DeepONet deep operator network. The corresponding plastic zone range and loose zone thickness are statistically analyzed by combining the Mohr-Coulomb yield criterion, and the surrounding rock state evolution sequence is constructed according to the excavation sequence. A 3D convolutional visual Transformer encoder is used to encode the spatial features of the three-dimensional physical field at each excavation length, obtain the corresponding high-dimensional feature vectors, and arrange them according to the excavation sequence. The excavation step sequence attention matrix is calculated to identify the key turning excavation lengths in the evolution of the surrounding rock state. The PCMCI causal discovery algorithm is run with the excavation step sequence attention matrix as a prior constraint. False causal edges caused by common dependent variables are eliminated through conditional independence tests. A directed acyclic causal graph containing key state variables and directed causal edges is output. The causal evolution path of the surrounding rock disaster along the excavation step sequence is quantitatively extracted from the directed acyclic causal graph.