Deep layered cataclastic rock large deformation analysis method and system based on multi-source fusion

Abstract

The application discloses a deep layered cataclastic rock large deformation analysis method and system based on multi-source fusion. By fusing geological survey and three-dimensional laser scanning data, the geometric characteristics of the structural plane are extracted to construct a discrete fracture network model, and the size of the rock mass representation unit (REV) is determined. Accordingly, a synthetic rock mass model with an elastic boundary is constructed, a three-dimensional cohesive block damage constitutive model is coupled, a linear degradation law of joint parameters with inelastic deformation is defined, and a step-by-step excavation model is established. The measured ground stress is applied, and the stress pullback correction algorithm based on the vertical partition of the plastic potential function is used for iterative calculation, the damage parameters are updated in real time, and the displacement field and damage distribution are output. Finally, according to the damage characteristics and displacement vector, the V-shaped damage zone of the top arch, the extrusion zone of the side wall and the rock block dislocation zone are accurately identified. The application solves the problems of difficulty in representing the randomness of the cataclastic rock structural plane and inaccuracy in simulating the dynamic degradation of parameters, and significantly improves the large deformation analysis accuracy.

Description

Technical Field

[0001] This invention relates to the field of geological disaster simulation and early warning technology, and in particular to a method and system for analyzing large deformations of deep layered fractured rocks based on multi-source fusion. Background Technology

[0002] Layered fractured rock tunnels in deep mines, transportation tunnels, and hydropower station caverns under high stress environments are prone to progressive fracturing and large deformation disasters after excavation and unloading due to the high strength of the rock blocks, well-developed joints, and low rock mass strength. This seriously threatens construction and operational safety. Current synthetic rock mass models provide an effective way to analyze the mechanical response of fractured rock masses by combining discrete fracture networks with numerical simulation. However, the coupled joint constitutive models are mostly based on ideal elastoplastic models, which cannot simultaneously characterize the stiffness degradation, strength degradation, and damage accumulation laws of joints under tensile, compressive-shear, and tensile-shear combined stress states. It is also difficult to reproduce the progressive damage process of complex stress paths on structural surfaces under high stress.

[0003] Given the unique mechanical properties of layered fractured hard rock, traditional constitutive models struggle to differentiate between various failure modes, such as tension, compression-shear, and tensile-shear, and cannot quantitatively output damage distribution, making it difficult to support precise support design based on simulation results. Furthermore, existing methods for simulating the excavation process often employ a single, complete excavation, failing to recreate the evolution of surrounding rock damage under multi-stage unloading and thus failing to accurately reflect the development patterns of large deformations controlled by structural planes.

[0004] Therefore, there is an urgent need to develop analytical methods and systems that can quantitatively characterize progressive damage to structural surfaces and distinguish failure modes, so as to provide a reliable basis for stability assessment and support optimization of deep fractured rock tunnels. Summary of the Invention

[0005] To address at least one of the aforementioned technical problems, this invention proposes a method and system for analyzing large deformations in deep layered fractured rocks based on multi-source fusion.

[0006] The first aspect of this invention provides a method for analyzing large deformations in deep layered fractured rocks based on multi-source fusion, comprising: Geological exploration data and three-dimensional laser scanning point cloud data of deep layered fractured rock tunnels are obtained. Based on the three-dimensional laser scanning point cloud data, the geometric features of rock mass structural surfaces are extracted to obtain a geometric feature dataset of structural surfaces. The geometric feature dataset of the structural surface is reconstructed based on the discrete fracture network algorithm to obtain a discrete fracture network model. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterization unit is determined. A numerical model of synthetic rock mass with elastic boundaries is constructed using the unit size of the rock mass characterization unit. A three-dimensional bonded block damage constitutive model is coupled with it. The linear deterioration law of joint normal stiffness, tangential stiffness, tensile strength, cohesion and dilatation angle with inelastic deformation is defined. The surrounding rock units are layered and grouped according to the excavation steps to form a step-by-step excavation simulation model. The measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and the excavation calculation is carried out step by step. In each calculation step, when the stress of the joint unit reaches the yield condition, stress pull-back correction is performed based on the vertical partition of the plastic potential function, the inelastic deformation and joint mechanical parameters are iteratively updated, and the surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Based on the distribution of damage types, damage degree, and displacement vectors, the V-shaped failure zone of the top arch, the extrusion zone of the sidewalls, and the fault zone of the layered rock blocks are identified, and the large deformation analysis results of the deep layered fractured rock are obtained.

[0007] In this scheme, the acquisition of geological exploration data and three-dimensional laser scanning point cloud data of deep layered fractured rock tunnels, and the extraction of geometric features of rock mass structural surfaces based on the three-dimensional laser scanning point cloud data to obtain a structural surface geometric feature dataset, specifically involves: Geological survey data of deep layered fractured rock tunnels are obtained, stratigraphic lithology information, joint statistics and geological structure distribution information of the target area are extracted, and three-dimensional laser scanning equipment is used to perform full-area scanning of the tunnel face and exposed rock surface to obtain point cloud data of the rock surface. The point cloud data is subjected to noise reduction filtering and coordinate registration to construct a three-dimensional digital model of the rock mass surface. The three-dimensional digital model of the rock mass surface is then fitted with a local surface to calculate the normal vector of each data point in the point cloud. Based on the abrupt change feature of normal vector change, the structural surface traces of the rock mass are identified. The continuous structural surface boundary is traced along the trace direction. The closed structural surface boundary is divided into regions to obtain planar fragments of individual structural surfaces. Spatial plane fitting is performed on each planar fragment to calculate the spatial dip and dip angle of the planar fragment. The spacing between adjacent structural surface traces in the same plane is statistically analyzed. Multiple groups of planar fragments are grouped according to clustering thresholds of dip and tilt angle. The average spacing of structural surfaces and the number of structural surfaces per unit area in each group are calculated to obtain the density parameters of each group of structural surfaces. The spatial geometric parameters of structural surfaces in all groups are summarized to construct a structural surface geometric feature dataset containing structural surface dip, tilt angle, spacing, density, and distribution range.

[0008] In this scheme, the geometric feature dataset of the structural surface is reconstructed using a discrete fracture network algorithm to obtain a discrete fracture network model. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterizing unit is determined. Specifically: Using the dip, tilt angle, spacing and density parameters in the geometric feature dataset of the structural surface as input variables, the Monte Carlo random sampling method is used to generate multiple sets of independent random structural surfaces in the three-dimensional spatial domain. The intersection operation between each set of random structural surfaces and the boundary of the spatial domain is performed to extract the trace distribution features of the structural surface on the boundary. Based on the trace distribution characteristics, a three-dimensional discrete fracture network model consistent with the orientation distribution of the structural surface in the target region is constructed. Based on the three-dimensional discrete fracture network model, multiple cubic rock mass models are set in order of increasing size. The discrete fracture network in each cubic rock mass model is spatially cut to obtain a set of fragmented rock blocks of the corresponding size. Apply uniaxial compression boundary conditions to each set of fractured rock blocks, calculate the equivalent compressive strength of the set of fractured rock blocks at each size, and plot the convergence curve of the equivalent compressive strength as the size of the cubic rock mass model changes. When the change in equivalent compressive strength during the increase of the size of the cubic rock mass model is lower than the set convergence threshold, the size is determined as the size of the rock mass characterization unit.

[0009] In this scheme, a synthetic rock mass numerical model with elastic boundaries is constructed using the dimensions of the rock mass characterizing unit. This model is coupled with a three-dimensional bonded block damage constitutive model. The linear degradation laws of joint normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle with inelastic deformation are defined. The surrounding rock units are layered and grouped according to the excavation sequence to form a step-by-step excavation simulation model. Specifically: The analysis range of the discrete fracture network is defined based on the size of the rock mass characterization unit. All structural surfaces in the discrete fracture network are retained within the analysis range. The structural surfaces are used as cutting boundaries to spatially divide the rock blocks, resulting in a discretely distributed set of fractured rock blocks. Elastic boundary blocks with a size larger than the analysis range are generated at the critical periphery of the analysis range. The mechanical parameters of the elastic boundary blocks are consistent with the mechanical parameters of the intact rock blocks to eliminate the influence of boundary effects on the calculation results. The tunnel space is cut out inside the fragmented rock block set according to the design cross-sectional shape and size of the tunnel. The surrounding rock units within a predetermined range around the tunnel are divided into multiple excavation groups according to the actual excavation advance and construction sequence. Each excavation group corresponds to an independent excavation stage. The constitutive model of the intact rock block in the set of fractured rock blocks is set as an elastic constitutive model, and the constitutive model of the structural surface unit is set as a three-dimensional bonded block damage constitutive model. The three-dimensional bonded block damage constitutive model decomposes the displacement vector of the structural surface into elastic deformation part, plastic deformation part and fracture deformation part, and merges the plastic deformation part and fracture deformation part into inelastic deformation amount under tension, tension-shear and compression-shear stress states. Once the structural surface reaches the yield condition, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle are defined as linear functions of the inelastic deformation. Under compression, the structural surface exhibits only linear elastic behavior and does not enter the yield stage. Under tension and tension-shear conditions, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle all decrease linearly with the increase of inelastic deformation. Under compression-shear conditions, only tensile strength, cohesion, and dilatation angle decrease linearly with the increase of inelastic deformation. The normal stiffness and tangential stiffness remain unchanged from their initial values. After constitutive coupling and element grouping of the model, a step-by-step excavation simulation model containing elastic boundaries, fractured rock blocks, structural surfaces, and layered excavation elements is obtained.

[0010] In this scheme, the measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and the excavation calculation is carried out step by step. In each calculation step, when the stress of the joint element reaches the yield condition, stress pullback correction is performed based on the vertical partition of the plastic potential function, and the inelastic deformation and joint mechanical parameters are iteratively updated. The surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Specifically: An initial geostress field measured on-site is applied to the step-by-step excavation simulation model to constrain the normal displacement of the model boundary. The empty constitutive structure of the corresponding surrounding rock unit is activated in the order of the excavation group to simulate the excavation unloading effect. In each calculation step, the maximum unbalanced force of the model tends to stabilize as the convergence criterion. Calculate the elastic test stress of the joint element, substitute the elastic test stress into the yield surface function for judgment, when the elastic test stress meets the yield condition, construct the plastic potential function under the current calculation step, the plastic potential function is the linear relationship between tangential stress and normal stress, draw a perpendicular line to the plastic potential function through the origin, and use the perpendicular line as the partition boundary to divide the stress region above the yield surface into the lower left region and the upper right region. For the elastic test stress of the joint element located in the lower left region, calculate the slope of the line connecting the elastic test stress and the origin, project the elastic test stress onto the yield surface along the direction of the line to obtain the corrected joint element stress, and extract the components of the direction of the line in the normal and tangential directions as inelastic strain components. For the elastic test stress of the joint element located in the upper right region, calculate the normal vector of the plastic potential function at the origin, project the elastic test stress onto the yield surface along the direction of the normal vector to obtain the corrected joint element stress, and calculate the inelastic strain components based on the components of the normal vector in the normal and tangential directions. Divide the inelastic strain components by the initial normal stiffness and the initial tangential stiffness to obtain the inelastic strain coefficients in the normal and tangential directions, respectively. Combine the corrected joint element stress to calculate the inelastic multiplier. Multiply the inelastic multiplier by the inelastic strain coefficients in the normal and tangential directions respectively and sum them to obtain the inelastic deformation increment of the current calculation step. Add the inelastic deformation amount of the previous calculation step to obtain the total inelastic deformation amount of the current calculation step. The normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle of the joint are updated according to the total inelastic deformation in accordance with the linear deterioration law. The damage degree of the bonded block is calculated based on the ratio of the updated normal stiffness to the initial normal stiffness. The displacement values ​​of each node of the surrounding rock are obtained by solving the overall equilibrium equation, and the displacement field of the surrounding rock is constructed. The damage type of the structural surface is determined according to the stress state and inelastic deformation component of the modified joint element. The damage type includes tensile damage, compression-shear damage and tension-shear damage. The surrounding rock displacement field, the distribution of structural surface damage types and the corresponding damage degree data are output simultaneously.

[0011] In this scheme, based on the distribution of damage types, damage degree, and displacement vectors, the V-shaped failure zone of the crown arch, the extrusion zone of the sidewalls, and the faulting zone of layered rock blocks are identified, resulting in the analysis results of large deformation of deep layered fractured rock. Specifically: The displacement vector of each node in the surrounding rock displacement field output by each calculation step of the step-by-step excavation simulation model is retrieved. The damage type and corresponding damage degree of the joint unit at the location of each node are associated. The arch area and the sidewall area are divided with the center line of the roadway as the reference. The damage type of each joint unit in the arch area is traversed. When the damage type of the joint unit is tensile damage and the damage degree reaches the maximum value, and the displacement vector of the node corresponding to the joint unit points to the free face of the roadway, the area where the joint unit is located is marked as the arch V-shaped failure zone. The damage type of each joint unit in the sidewall region is traversed. When the damage type of the joint unit is compressive-shear damage and the damage degree reaches the maximum value, and the displacement vector of the corresponding node of the joint unit points to the free surface inside the roadway, the area where the joint unit is located is marked as the sidewall extrusion zone. Paired tracking of joint units on both sides of the layered structural plane. When both paired joint units exhibit tensile-shear damage and the damage degree continues to increase with the time step, and the displacement vectors of the corresponding nodes of the two joint units have relative shear displacement along the direction of the structural plane, the area where the structural plane is located is marked as the layered rock block displacement zone. By summarizing the spatial distribution range, corresponding damage extreme values, and displacement vector amplitudes of all marked V-shaped failure zones of the top arch, extrusion zones of the sidewalls, and fault zones of layered rock blocks, a mapping relationship between large deformation modes and spatial locations is constructed, and the large deformation analysis results of deep layered fractured rocks are obtained.

[0012] A second aspect of the present invention provides a system for analyzing large deformations of deep layered fractured rocks based on multi-source fusion. The system includes a memory and a processor. The memory includes a program for analyzing large deformations of deep layered fractured rocks based on multi-source fusion. When the program for analyzing large deformations of deep layered fractured rocks based on multi-source fusion is executed by the processor, it implements the steps of the method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion as described in any of the preceding claims.

[0013] This invention discloses a method and system for large deformation analysis of deep layered fractured rock based on multi-source fusion. By fusing geological survey and 3D laser scanning data, geometric features of structural surfaces are extracted to construct a discrete fracture network model, and the dimensions of the rock mass characterization unit (REV) are determined. Based on this, a synthetic rock mass model with elastic boundaries is constructed, coupled with a 3D bonded block damage constitutive model, defining the linear degradation law of joint parameters with inelastic deformation, and establishing a step-by-step excavation model. Measured in-situ stress is applied, and a stress pullback correction algorithm based on vertical partitioning of plastic potential function is used for iterative calculation, updating damage parameters in real time, and outputting the displacement field and damage distribution. Finally, based on damage characteristics and displacement vectors, the V-shaped failure zone of the crown arch, the extrusion zone of the sidewalls, and the rock block displacement zone are accurately identified. This invention solves the problems of difficulty in characterizing the randomness of fractured rock structural surfaces and inaccurate simulation of dynamic parameter degradation, significantly improving the accuracy of large deformation analysis. Attached Figure Description

[0014] Figure 1 A flowchart of a method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to the present invention is shown; Figure 2 This is a schematic diagram of stress pullback correction under different joint element stress states in the three-dimensional bonded block damage constitutive model of the present invention; Figure 3 A schematic diagram of structural surface damage types under different joint element stress states in a three-dimensional bonded block damage constitutive model; Wherein, 1-elastic test stress of joint element, 2-stress state of joint element after stress pullback correction, 3-yield surface under the current calculation step, 4-plastic potential function under the current calculation step, 5-stress state of joint under tensile damage state, 6-stress state of joint under tensile-shear damage state, 7-stress state of joint under compressive-shear state, 8-stress state of joint under tensile-shear action, 9-stress state of joint under compressive-shear action; Figure 4 A block diagram of a deep layered fractured rock large deformation analysis system based on multi-source fusion is shown. Detailed Implementation

[0015] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0016] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0017] Figure 1 The flowchart of a method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to the present invention is shown.

[0018] like Figure 1 As shown, the first aspect of the present invention provides a method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion, comprising: Geological exploration data and three-dimensional laser scanning point cloud data of deep layered fractured rock tunnels are obtained. Based on the three-dimensional laser scanning point cloud data, the geometric features of rock mass structural surfaces are extracted to obtain a geometric feature dataset of structural surfaces. The geometric feature dataset of the structural surface is reconstructed based on the discrete fracture network algorithm to obtain a discrete fracture network model. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterization unit is determined. A numerical model of synthetic rock mass with elastic boundaries is constructed using the unit size of the rock mass characterization unit. A three-dimensional bonded block damage constitutive model is coupled with it. The linear deterioration law of joint normal stiffness, tangential stiffness, tensile strength, cohesion and dilatation angle with inelastic deformation is defined. The surrounding rock units are layered and grouped according to the excavation steps to form a step-by-step excavation simulation model. The measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and the excavation calculation is carried out step by step. In each calculation step, when the stress of the joint unit reaches the yield condition, stress pull-back correction is performed based on the vertical partition of the plastic potential function, the inelastic deformation and joint mechanical parameters are iteratively updated, and the surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Based on the distribution of damage types, damage degree, and displacement vectors, the V-shaped failure zone of the top arch, the extrusion zone of the sidewalls, and the fault zone of the layered rock blocks are identified, and the large deformation analysis results of the deep layered fractured rock are obtained.

[0019] It should be noted that by fusing multi-source data, the natural structural features of deep layered fractured rocks are accurately characterized. Three-dimensional laser scanning point clouds are used to identify the geometric parameters of structural surfaces and construct a discrete fracture network, solving the problem that traditional exploration methods struggle to reconstruct the spatial distribution of complex joints and ensuring the accuracy of the model's geological prototype. The size of the characterizing unit is determined by the convergence of equivalent compressive strength with rock mass size, eliminating the influence of size effects and enabling the subsequent numerical model to represent the real mechanical behavior of the rock mass. A synthetic rock mass model with elastic boundaries is constructed and coupled with a three-dimensional bonded block damage constitutive model, overcoming the limitation of traditional ideal elastoplastic models in distinguishing between different failure mechanisms such as tension, compression-shear, and tension-shear. By defining the linear degradation law of joint mechanical parameters with inelastic deformation, the entire process of structural surfaces from crack initiation and damage to complete failure is realized. The method employs quantitative characterization; it simulates the unloading process in layers and groups according to the excavation steps, restoring the stress path evolution characteristics caused by step excavation and avoiding the distortion of the damage process caused by one-time excavation; based on the plastic potential function, it performs stress pullback correction by vertical partitioning, accurately capturing the yield response of joint units under different stress states, realizing dynamic iterative updates of inelastic deformation and damage degree, and ensuring that the output of surrounding rock displacement field and structural surface damage type has both physical mechanism rationality and calculation accuracy; finally, through the joint identification of damage type, damage degree and displacement vector, it accurately distinguishes three typical large deformation modes: the V-shaped failure zone of the top arch, the extrusion zone of the sidewall, and the fault zone of layered rock blocks, revealing the evolution mechanism of structural surface-controlled large deformation, and providing a reliable quantitative basis for the stability assessment and support design of deep fractured rock tunnels. This method achieves refined analysis of the entire process of large deformation in deep layered fractured rock through full-chain technological innovation, from accurate reconstruction of geological structure, equivalent rock mass mechanical behavior, simulation of damage evolution process to intelligent identification of large deformation modes.

[0020] According to an embodiment of the present invention, the step of acquiring geological exploration data and three-dimensional laser scanning point cloud data of deep layered fractured rock tunnels, and extracting geometric features of rock mass structural surfaces based on the three-dimensional laser scanning point cloud data to obtain a structural surface geometric feature dataset, specifically includes: Geological survey data of deep layered fractured rock tunnels are obtained, stratigraphic lithology information, joint statistics and geological structure distribution information of the target area are extracted, and three-dimensional laser scanning equipment is used to perform full-area scanning of the tunnel face and exposed rock surface to obtain point cloud data of the rock surface. The point cloud data is subjected to noise reduction filtering and coordinate registration to construct a three-dimensional digital model of the rock mass surface. The three-dimensional digital model of the rock mass surface is then fitted with a local surface to calculate the normal vector of each data point in the point cloud. Based on the abrupt change feature of normal vector change, the structural surface traces of the rock mass are identified. The continuous structural surface boundary is traced along the trace direction. The closed structural surface boundary is divided into regions to obtain planar fragments of individual structural surfaces. Spatial plane fitting is performed on each planar fragment to calculate the spatial dip and dip angle of the planar fragment. The spacing between adjacent structural surface traces in the same plane is statistically analyzed. Multiple groups of planar fragments are grouped according to clustering thresholds of dip and tilt angle. The average spacing of structural surfaces and the number of structural surfaces per unit area in each group are calculated to obtain the density parameters of each group of structural surfaces. The spatial geometric parameters of structural surfaces in all groups are summarized to construct a structural surface geometric feature dataset containing structural surface dip, tilt angle, spacing, density, and distribution range.

[0021] It should be noted that the stratigraphic lithological information includes lithological type, single layer thickness, weathering degree, mineral composition and physical and mechanical parameters of rock blocks; joint statistics include structural surface type, infill material properties, roughness coefficient, opening degree and penetration rate; and geological structure distribution information includes fault attitude, fold morphology, fracture zone width and principal direction and magnitude of tectonic stress field.

[0022] According to an embodiment of the present invention, the discrete fracture network model is obtained by reconstructing the geometric feature dataset of the structural surface based on the discrete fracture network algorithm. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterizing unit is determined. Specifically: Using the dip, tilt angle, spacing and density parameters in the geometric feature dataset of the structural surface as input variables, the Monte Carlo random sampling method is used to generate multiple sets of independent random structural surfaces in the three-dimensional spatial domain. The intersection operation between each set of random structural surfaces and the boundary of the spatial domain is performed to extract the trace distribution features of the structural surface on the boundary. Based on the trace distribution characteristics, a three-dimensional discrete fracture network model consistent with the orientation distribution of the structural surface in the target region is constructed. Based on the three-dimensional discrete fracture network model, multiple cubic rock mass models are set in order of increasing size. The discrete fracture network in each cubic rock mass model is spatially cut to obtain a set of fragmented rock blocks of the corresponding size. Apply uniaxial compression boundary conditions to each set of fractured rock blocks, calculate the equivalent compressive strength of the set of fractured rock blocks at each size, and plot the convergence curve of the equivalent compressive strength as the size of the cubic rock mass model changes. When the change in equivalent compressive strength during the increase of the size of the cubic rock mass model is lower than the set convergence threshold, the size is determined as the size of the rock mass characterization unit.

[0023] It should be noted that the discrete fracture network model is constructed to overcome the limitations of traditional continuous medium models that treat rock mass as a homogeneous medium. By introducing the geometric characteristic parameters of structural surfaces measured in the field, the spatial distribution, connectivity, and combination characteristics of structural surfaces in layered fractured rock are accurately reproduced in the three-dimensional spatial domain, thereby truly reflecting the control effect of structural surfaces on the mechanical behavior of rock mass. Determining the size of the rock mass characterization unit is to solve the size effect problem of rock mass mechanical parameters. Since the equivalent compressive strength of fractured rock mass exhibits significant nonlinear characteristics as the research scale changes, the rock mass mechanical parameters will tend to stabilize only when the model size reaches the characterization unit scale. At this time, the model can contain a sufficient number of structural surfaces to reflect the statistical homogeneity of the rock mass, and can avoid computational redundancy caused by an excessively large model. This provides a basic analysis unit with both physical realism and computational feasibility for subsequent numerical simulation of synthetic rock mass, ensuring the reliability and engineering applicability of the large deformation analysis results. The uniaxial compression boundary condition refers to the combination of loading and constraint that applies a vertically downward axial pressure to the upper surface of the cubic rock mass model, completely fixes the lower surface, and removes the normal constraints on the other four sides to allow free deformation. Specifically, it includes axial displacement or load control methods, lateral zero normal stress boundaries, and rigid loading plate contact settings on the upper and lower end faces.

[0024] According to an embodiment of the present invention, the method of constructing a synthetic rock mass numerical model with elastic boundaries based on the size of the rock mass characterization unit, coupling a three-dimensional bonded block damage constitutive model, defining the linear degradation law of joint normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle with inelastic deformation, and grouping the surrounding rock units into layers according to the excavation sequence to form a step-by-step excavation simulation model, specifically: The analysis range of the discrete fracture network is defined based on the size of the rock mass characterization unit. All structural surfaces in the discrete fracture network are retained within the analysis range. The structural surfaces are used as cutting boundaries to spatially divide the rock blocks, resulting in a discretely distributed set of fractured rock blocks. Elastic boundary blocks with a size larger than the analysis range are generated at the critical periphery of the analysis range. The mechanical parameters of the elastic boundary blocks are consistent with the mechanical parameters of the intact rock blocks to eliminate the influence of boundary effects on the calculation results. The tunnel space is cut out inside the fragmented rock block set according to the design cross-sectional shape and size of the tunnel. The surrounding rock units within a predetermined range around the tunnel are divided into multiple excavation groups according to the actual excavation advance and construction sequence. Each excavation group corresponds to an independent excavation stage. The constitutive model of the intact rock block in the set of fractured rock blocks is set as an elastic constitutive model, and the constitutive model of the structural surface unit is set as a three-dimensional bonded block damage constitutive model. The three-dimensional bonded block damage constitutive model decomposes the displacement vector of the structural surface into elastic deformation part, plastic deformation part and fracture deformation part, and merges the plastic deformation part and fracture deformation part into inelastic deformation amount under tension, tension-shear and compression-shear stress states. Once the structural surface reaches the yield condition, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle are defined as linear functions of the inelastic deformation. Under compression, the structural surface exhibits only linear elastic behavior and does not enter the yield stage. Under tension and tension-shear conditions, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle all decrease linearly with the increase of inelastic deformation. Under compression-shear conditions, only tensile strength, cohesion, and dilatation angle decrease linearly with the increase of inelastic deformation. The normal stiffness and tangential stiffness remain unchanged from their initial values. After constitutive coupling and element grouping of the model, a step-by-step excavation simulation model containing elastic boundaries, fractured rock blocks, structural surfaces, and layered excavation elements is obtained.

[0025] It should be noted that the core features of the three-dimensional bonded block damage constitutive model are as follows: (1) Basic definition of the three-dimensional bonded block damage constitutive model This constitutive model decomposes the displacement vectors of the structural surfaces into elastic deformation, plastic deformation, and fracture deformation components, namely: (1) In the formula, For elastic deformation, For plastic deformation, This is called fracture deformation. Plastic deformation and fracture deformation are collectively referred to as inelastic deformation. .

[0026] (2) Mechanical behavior under different stress states The mechanical behavior of this constitutive model under different stress states is as follows: ① Under compression: the joint exhibits linear elasticity under compression and will not reach the yield state, i.e., no damage occurs; ② Under tension: under tension, the joint is linear elastic before yielding; after yielding, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle deteriorate linearly with the increase of inelastic deformation; ③ Under tension-shear: similar to the tensile state, after yielding, each parameter deteriorates linearly with inelastic deformation; ④ Under compression-shear: under compression-shear, after yielding, the tensile strength, cohesion, and dilatation angle deteriorate linearly with inelastic deformation, but the normal stiffness and tangential stiffness do not deteriorate.

[0027] Among them, the yield state refers to the critical state when the combination of normal stress and tangential stress of the joint element reaches the predefined strength envelope, i.e., the yield surface function, during the stress process. The joint element is the computational unit formed by discretizing the structural surface in the discrete fracture network in the synthetic rock mass numerical model, which is used to simulate the mechanical behavior of discontinuous surfaces inside the rock mass.

[0028] To achieve the aforementioned mechanical behavior, in the simulation and iterative calculation of the three-dimensional bonded block damage constitutive model, when the stress state of the joint element has not reached the yield state, the normal and tangential stress increments on the joint surface and the strain increments in the corresponding directions exhibit linear elastic changes, as shown in Equations 2 and 3: (2) (3) In the formula, , The normal and tangential stiffness are given in the current calculation step. , This represents the normal strain increment of the structural surface in the current calculation step. , This represents the increment of normal stress on the structural surface in the current calculation step.

[0029] Under the above conditions, when the joint stress state ( , When the yield condition is reached through iterative processes, the joint stress state calculated according to the linear elastic variation law will be pulled back to the yield surface, as shown in Equations 4 to 5.

[0030] (4) (5) In the formula, It is an inelastic multiplier. , The normal and tangential stiffness under initial conditions. , For the inelastic strain components in the normal and tangential directions.

[0031] (3) Parameter degradation law Under the influence of inelastic deformation, joints in tensile, tensile-shear, and compressive-shear stress states will eventually yield, leading to gradual damage to the structural surface. To illustrate the damage process of the structural surface, the tensile strength of the fracture surface is linearly degraded according to the Type I fracture energy, while the cohesion and dilatation angle are linearly degraded according to the Type II fracture energy. (6) (7) (8) (9) (10) (11) In the formula, This is the amount of inelastic deformation. , and Tensile strength Cohesion and shear expansion angle initial value, , and These represent the tensile strength, cohesion, and inelastic deformation values ​​when the shear dilatation angle is reduced to 0, respectively. It is a type I fracture energy. This is an approximate value for the type II fracture energy.

[0032] (4) Stiffness deterioration after yielding In the yield state, the normal and tangential stiffness between the bonded blocks are degraded. At this point, the stress state lies on the yield surface, and the normal stress is the degraded tensile strength. The tangential stress is calculated using the yield surface function, and the degraded normal stiffness is... With tangential stiffness Calculate according to Equations 12 and 13: (12) (13) In the formula, and For the initial normal stiffness and initial tangential stiffness, The plastic deformation portion of inelastic deformation (i.e.) ), This represents the joint shear strength under the current inelastic deformation state.

[0033] In the three-dimensional bonded block damage constitutive model, the joints are linearly elastic and do not yield under compression; under tension, tension-shear, and compression-shear conditions, they are linearly elastic before yielding, and after yielding, the normal stiffness, tangential stiffness, tensile strength, cohesion, and shear dilatation angle deteriorate linearly with the increase of inelastic deformation.

[0034] According to an embodiment of the present invention, the measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and the excavation calculation is carried out step by step. In each calculation step, when the stress of the joint element reaches the yield condition, stress pullback correction is performed based on the vertical partition of the plastic potential function, the inelastic deformation and joint mechanical parameters are iteratively updated, and the surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Specifically: An initial geostress field measured on-site is applied to the step-by-step excavation simulation model to constrain the normal displacement of the model boundary. The empty constitutive structure of the corresponding surrounding rock unit is activated in the order of the excavation group to simulate the excavation unloading effect. In each calculation step, the maximum unbalanced force of the model tends to stabilize as the convergence criterion. Calculate the elastic test stress of the joint element, substitute the elastic test stress into the yield surface function for judgment, when the elastic test stress meets the yield condition, construct the plastic potential function under the current calculation step, the plastic potential function is the linear relationship between tangential stress and normal stress, draw a perpendicular line to the plastic potential function through the origin, and use the perpendicular line as the partition boundary to divide the stress region above the yield surface into the lower left region and the upper right region. For the elastic test stress of the joint element located in the lower left region, calculate the slope of the line connecting the elastic test stress and the origin, project the elastic test stress onto the yield surface along the direction of the line to obtain the corrected joint element stress, and extract the components of the direction of the line in the normal and tangential directions as inelastic strain components. For the elastic test stress of the joint element located in the upper right region, calculate the normal vector of the plastic potential function at the origin, project the elastic test stress onto the yield surface along the direction of the normal vector to obtain the corrected joint element stress, and calculate the inelastic strain components based on the components of the normal vector in the normal and tangential directions. Divide the inelastic strain components by the initial normal stiffness and the initial tangential stiffness to obtain the inelastic strain coefficients in the normal and tangential directions, respectively. Combine the corrected joint element stress to calculate the inelastic multiplier. Multiply the inelastic multiplier by the inelastic strain coefficients in the normal and tangential directions respectively and sum them to obtain the inelastic deformation increment of the current calculation step. Add the inelastic deformation amount of the previous calculation step to obtain the total inelastic deformation amount of the current calculation step. The normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle of the joint are updated according to the total inelastic deformation in accordance with the linear deterioration law. The damage degree of the bonded block is calculated based on the ratio of the updated normal stiffness to the initial normal stiffness. The displacement values ​​of each node of the surrounding rock are obtained by solving the overall equilibrium equation, and the displacement field of the surrounding rock is constructed. The damage type of the structural surface is determined according to the stress state and inelastic deformation component of the modified joint element. The damage type includes tensile damage, compression-shear damage and tension-shear damage. The surrounding rock displacement field, the distribution of structural surface damage types and the corresponding damage degree data are output simultaneously.

[0035] It should be noted that the initial geostress field obtained by applying measured or inverted data to the step-by-step excavation simulation model includes three normal stress components. , , and three shear stress components , , The normal displacement on the model boundary is constrained, ensuring that the velocity vector of the elements on the boundary remains zero throughout the simulation. Joint mechanical parameters include the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle at the initial and current calculation steps, as well as intermediate values ​​during the linear degradation to zero with inelastic deformation and the degree of damage to the bonded block calculated based on the stiffness degradation ratio.

[0036] The simulation process is carried out step by step according to the excavation stages. In each excavation stage, the corresponding group of roadway elements are given empty constitutive models and numerical calculations are performed to simulate the unloading effect of the roadway segment excavation process. The number of calculation steps in each stage is converged when the maximum unbalanced force of the model tends to stabilize.

[0037] In each computational step, the following constitutive update and computation process is executed: Figure 2 This is a schematic diagram of stress pullback correction under different joint element stress states in the three-dimensional bonded block damage constitutive model of the present invention.

[0038] (1) Stress pullback correction When the stress state of the joint element reaches the yield condition, a partitioned stress pull-back algorithm is used to correct the elastic test stress to the yield surface. During this process, for the inelastic strain components... and The calculation process is determined based on the stress state of the joint element at yield. The plastic potential and yield surface function are defined by equations 15 and 16. A perpendicular line OP is drawn from the origin to the plastic potential function, which divides the region above the yield surface into two parts: the lower left and the upper right. For the stress state of the joint element on the lower left side of OP, it is pulled back to the yield surface through the direction of the line connecting it to the origin; while for the stress state of the joint element on the upper right side of OP, it is pulled back through a direction perpendicular to the plastic potential function. Thus, the corresponding formulas for calculating the inelastic strain components are shown in equations 16 and 17: (14) (15) (16) (17) (2) Calculation of inelastic deformation and damage After obtaining the inelastic strain components, substituting them into Equations 4 and 5 allows for the further solution of the inelastic multipliers. Thus, the inelastic deformation in the current N+1th calculation step is calculated according to Equation 18.

[0039] (18) The degree of damage to the bonded mass is defined based on the degree of degradation of the normal stiffness during the calculation process. Used to indicate the degree of damage between blocks. The inelastic deformation is Tensile strength at that time This refers to the plastic deformation portion within the inelastic deformation. For the initial normal stiffness, For the normal component of elastic deformation.

[0040] . The range of change is 0-1, and the larger the value, the greater the degree of damage. The calculation is performed according to Equation 19: (19) (3) Output of surrounding rock displacement field and structural surface damage type Figure 3 Schematic diagram of structural surface damage types under different joint element stress states in a three-dimensional bonded block damage constitutive model. After completing the constitutive update in this computational step, the displacement values ​​of each node in the surrounding rock are obtained by solving the global equilibrium equations, forming the surrounding rock displacement field. Simultaneously, based on the stress state and inelastic deformation components of each joint element, the damage type of the structural surface is determined. ① Tensile damage: The stress state of the joint element is The corresponding joint element yields under tension or tension-shear conditions; ② Compression-shear damage: The stress state of the joint element is , The corresponding joint element yields under compressive shear conditions; ③ Tensile-shear damage: The stress state of the joint element is , The corresponding joint element yields under tension and shear conditions.

[0041] The overall equilibrium equation is a matrix governing equation based on the principle of virtual work, characterizing the force balance of all nodes in the model during step-by-step excavation simulation. Its expression is K⋅U=R, where K is the overall stiffness matrix including the elastic stiffness of the intact rock block and the tangential stiffness of the joint element after stress pullback correction (deteriorating with inelastic deformation under tension / tension-shear conditions, maintaining the initial value of normal / tangential stiffness under compression-shear conditions, and linear elasticity under compression conditions), U is the displacement column vector of each node of the surrounding rock to be solved, and R is the external load column vector of the current calculation step—composed of the equivalent node force obtained by converting the initial ground stress measured on site and subtracting the load released by the corresponding excavation group element. During the equation solution process, the normal displacement constraint condition of the model boundary is satisfied at the same time. Through iteration, the maximum unbalanced force of the model tends to stabilize, and finally the displacement field and updated stress field of the surrounding rock under the current calculation step are obtained.

[0042] According to an embodiment of the present invention, the step of identifying the V-shaped failure zone of the crown arch, the extrusion zone of the sidewall, and the faulting zone of the layered rock blocks based on the damage type distribution, damage degree, and displacement vector, and obtaining the large deformation analysis results of the deep layered fractured rock, specifically includes: The displacement vector of each node in the surrounding rock displacement field output by each calculation step of the step-by-step excavation simulation model is retrieved. The damage type and corresponding damage degree of the joint unit at the location of each node are associated. The arch area and the sidewall area are divided with the center line of the roadway as the reference. The damage type of each joint unit in the arch area is traversed. When the damage type of the joint unit is tensile damage and the damage degree reaches the maximum value, and the displacement vector of the node corresponding to the joint unit points to the free face of the roadway, the area where the joint unit is located is marked as the arch V-shaped failure zone. The damage type of each joint unit in the sidewall region is traversed. When the damage type of the joint unit is compressive-shear damage and the damage degree reaches the maximum value, and the displacement vector of the corresponding node of the joint unit points to the free surface inside the roadway, the area where the joint unit is located is marked as the sidewall extrusion zone. Paired tracking of joint units on both sides of the layered structural plane. When both paired joint units exhibit tensile-shear damage and the damage degree continues to increase with the time step, and the displacement vectors of the corresponding nodes of the two joint units have relative shear displacement along the direction of the structural plane, the area where the structural plane is located is marked as the layered rock block displacement zone. By summarizing the spatial distribution range, corresponding damage extreme values, and displacement vector amplitudes of all marked V-shaped failure zones of the top arch, extrusion zones of the sidewalls, and fault zones of layered rock blocks, a mapping relationship between large deformation modes and spatial locations is constructed, and the large deformation analysis results of deep layered fractured rocks are obtained.

[0043] It should be noted that, based on the damage distribution and displacement values ​​of joint elements in each region and the damage type of the structural surface, the large deformation mode is determined according to the following criteria: (1) Identification of V-shaped failure zone of the top arch When tensile damage occurs in the crown arch area and the damage degree reaches 1.0, and the displacement vector of the surrounding rock in this area points towards the free face, it is identified as a V-shaped failure zone of the crown arch. The formation mechanism of this type of failure zone is as follows: after the tunnel is excavated, stress redistribution forms tangential compressive stress around the tunnel outline. The structural surface at the crown arch location is subjected to combined tensile and shear forces, and the tensile strength continues to deteriorate. Ultimately, this leads to the rock blocks falling or sliding towards the free face under the action of gravity, forming a V-shaped collapse zone.

[0044] (2) Identification of the extrusion zone of the sidewall When compressive-shear damage occurs in the sidewall area, with a damage degree reaching 1.0, and the displacement vector points towards the interior of the roadway (free face), it is identified as a sidewall extrusion zone. The formation mechanism of this type of failure zone is as follows: the layered structural plane is subjected to tangential compressive stress, the friction between the structural planes is insufficient to resist the shear force, resulting in shear slip, the cohesion and tensile strength of the structural planes continue to deteriorate, and the thin layered rock mass is continuously extruded into the roadway under the action of compressive stress.

[0045] (3) Identification of fault zones in layered rock blocks When tensile and shear damage occurs on both sides of a layered structural plane, the degree of damage continues to increase, and relative shear displacement occurs between the blocks on both sides, it is identified as a layered rock block displacement zone. The formation mechanism of this type of failure zone is as follows: under the combined action of tensile and shear stress, the structural plane simultaneously undergoes normal opening and tangential slip, the shear resistance of the structural plane continuously decreases, and the rock blocks detach and undergo displacement deformation along the structural plane.

[0046] According to an embodiment of the present invention, it further includes: The three-dimensional laser scanning point cloud data is preprocessed. The signal attenuation area is identified based on the echo intensity distribution characteristics of the lower part of the face in the point cloud. The elevation abrupt change characteristics of the point cloud elevation field are combined to lock the elevation anomaly accumulation zone. The overlapping area of ​​the signal attenuation area and the elevation anomaly accumulation zone is marked as the rock powder covered area. The boundary contour of the rock powder covered area and the elevation mean of the internal point cloud are extracted. The thickness of the rock powder accumulation is calculated based on the average elevation of the rock powder-covered area and the elevation difference between the surrounding exposed bedrock. The data in the geometric feature dataset of the structural surface located within the boundary contour of the rock powder-covered area are removed to obtain the effective geometric feature dataset of the structural surface after removing the interfering data. Joint statistics are retrieved from geological exploration data. Based on the occurrence information of the dominant structural planes extending below the rock powder-covered area, the dip, dip angle, and density parameters of the structural planes to be completed are determined. The effective structural plane geometric feature dataset is reconstructed using a discrete fracture network algorithm. During the reconstruction process, the dip, dip angle and density parameters of the structural surface to be completed are used as constraints. Virtual structural surface interpolation is performed within the boundary contour of the rock powder covered area so that the virtual structural surface generated by interpolation is smoothly connected with the visible structural surface at the edge of the boundary contour in terms of attitude and position, resulting in an enhanced discrete fracture network model containing hidden structural surfaces. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the enhanced discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterizing unit is determined. The assumption of no structure surface at the bottom is corrected by the enhanced discrete fracture network model.

[0047] According to an embodiment of the present invention, the network reconstruction of the effective structural surface geometric feature dataset based on the discrete fracture network algorithm specifically includes: Using the inclination, tilt angle, and average spacing of the dominant structural surface as statistical distribution inputs, a stochastic probability model is established for the bottom concealed region. The spatial domain of random sampling is defined according to the boundary range of the bottom concealed region, and the coordinates of the center point of the stochastic structural surface are generated within the spatial domain. Based on the normal vector direction of the dominant structural surface group, determine the spatial rotation matrix of the random structural surface, project the random structural surface onto the top surface of the bottom hidden area, and calculate the orientation angle between the random structural surface and the top visible structural surface. When the included angle of the orientation is less than a preset tolerance threshold, the random structural surface is retained as a virtual structural surface fragment, and the virtual structural surface fragment is topologically connected to the top visible structural surface to construct a continuous structural surface trace that runs through the rock powder overburden. The continuous structural surface trace is used to correct the fracture connectivity rate of the bottom region in the enhanced discrete fracture network model.

[0048] According to an embodiment of the present invention, the step of identifying the signal attenuation zone based on the echo intensity distribution characteristics of the lower region of the face in the point cloud, and locking the elevation anomaly accumulation zone by combining the gradient abrupt change characteristics of the point cloud elevation field, specifically includes: The three-dimensional laser scanning point cloud data is processed into a grid, and the average echo intensity of the point cloud in each grid cell is calculated. Areas with average echo intensity lower than the strength threshold of the surrounding intact rock wall and continuous distribution are identified as signal attenuation areas. Simultaneously, elevation difference calculation is performed on the grid cells, and regions with elevation change rates exceeding preset steep change thresholds are extracted as elevation anomaly accumulation zones. The signal attenuation zone and the elevation anomaly accumulation zone are spatially superimposed to retain the intersection region that simultaneously satisfies the conditions of signal attenuation and abrupt elevation change. The boundary contour of the rock powder covered area is determined based on the spatial coordinate extreme values ​​of the point cloud within the intersection region. The rock powder accumulation thickness is calculated based on the elevation difference between the highest and lowest points within the intersection region. The boundary contour and rock powder accumulation thickness of the rock powder covered area are used as constraints to limit the spatial range of subsequent virtual structural surface interpolation.

[0049] It should be noted that after blasting and excavation in deep tunnels, a thick layer of rock powder and gravel often accumulates at the bottom of the working face. This area is precisely the key exposure location of the gently dipping dominant structural plane at the bottom of the layered fractured rock. Conventional 3D laser scanning is affected by the physical coverage of rock powder and the strong scattering and attenuation of laser light in loose media, making it impossible to obtain information on the true structural plane at the bottom. This leads to the bottom area being incorrectly filled as a complete rock mass without structural planes when constructing a discrete fracture network model. Consequently, when calculating equivalent mechanical parameters, the size of the determined rock mass characterization unit is too small due to the neglect of potential slip surfaces, making it difficult to capture the slip and uplift characteristics of the bottom rock blocks, resulting in blind spots in the identification of the sidewall extrusion zone and the layered rock block displacement zone. This invention uses gridding processing of 3D laser scanning point cloud data, combined with the point cloud echo intensity distribution characteristics and the elevation field gradient abrupt change characteristics, to accurately identify signal attenuation zones and elevation anomaly accumulations. The intersection of accumulation zones is used to accurately define the boundary contour of the rock powder-covered area and calculate the accumulation thickness. Then, the distorted pseudo-structural surface data in this area are eliminated to obtain an effective dataset that reflects the real rock mass structure. Based on this, the attitude information of dominant structural surfaces in geological exploration data is used to perform constrained virtual structural surface interpolation in the rock powder-covered area. This makes the generated virtual structural surfaces smoothly connected with the surrounding visible structural surfaces in terms of attitude and spatial location, constructing a continuous structural surface trace that runs through the rock powder cover layer. This significantly improves the fracture connectivity rate in the bottom region of the enhanced discrete fracture network model, effectively corrects the bias of the no-structure-surface assumption caused by the lack of bottom data in traditional methods, and ensures the accuracy of the rock mass characterization unit size determination. This provides a reliable model basis for accurately reproducing the bottom plate bulge and sidewall compression large deformation in the subsequent step-by-step excavation simulation.

[0050] Figure 4 A block diagram of a deep layered fractured rock large deformation analysis system based on multi-source fusion is shown.

[0051] A second aspect of the present invention also provides a system for analyzing large deformations of deep layered fractured rocks based on multi-source fusion. The system includes: a memory 401, a processor 402, and a communication interface 403. The memory includes a program for analyzing large deformations of deep layered fractured rocks based on multi-source fusion. The communication interface is used for data connection and communication between the memory and the processor. When the program for analyzing large deformations of deep layered fractured rocks based on multi-source fusion is executed by the processor, it implements the steps of the method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion as described in any of the above claims.

[0052] This invention discloses a method and system for large deformation analysis of deep layered fractured rock based on multi-source fusion. By fusing geological survey and 3D laser scanning data, geometric features of structural surfaces are extracted to construct a discrete fracture network model, and the dimensions of the rock mass characterization unit (REV) are determined. Based on this, a synthetic rock mass model with elastic boundaries is constructed, coupled with a 3D bonded block damage constitutive model, defining the linear degradation law of joint parameters with inelastic deformation, and establishing a step-by-step excavation model. Measured in-situ stress is applied, and a stress pullback correction algorithm based on vertical partitioning of plastic potential function is used for iterative calculation, updating damage parameters in real time, and outputting the displacement field and damage distribution. Finally, based on damage characteristics and displacement vectors, the V-shaped failure zone of the crown arch, the extrusion zone of the sidewalls, and the rock block displacement zone are accurately identified. This invention solves the problems of difficulty in characterizing the randomness of fractured rock structural surfaces and inaccurate simulation of dynamic parameter degradation, significantly improving the accuracy of large deformation analysis.

[0053] Those skilled in the art will understand that all or part of the steps of the above method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it performs the steps of the above method embodiments. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0054] Alternatively, if the integrated units of this invention are implemented as software functional modules and sold or used as independent products, they can also be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the embodiments of this invention, or the parts that contribute to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as mobile storage devices, ROM, RAM, magnetic disks, or optical disks.

[0055] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion, characterized in that, Includes the following steps: Geological exploration data and three-dimensional laser scanning point cloud data of deep layered fractured rock tunnels are obtained. Based on the three-dimensional laser scanning point cloud data, the geometric features of rock mass structural surfaces are extracted to obtain a geometric feature dataset of structural surfaces. The geometric feature dataset of the structural surface is reconstructed based on the discrete fracture network algorithm to obtain a discrete fracture network model. The equivalent mechanical parameters of rock masses of different sizes are calculated based on the discrete fracture network model. When the equivalent mechanical parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterization unit is determined. A numerical model of synthetic rock mass with elastic boundaries is constructed using the unit size of the rock mass characterization unit. A three-dimensional bonded block damage constitutive model is coupled with it. The linear deterioration law of joint normal stiffness, tangential stiffness, tensile strength, cohesion and dilatation angle with inelastic deformation is defined. The surrounding rock units are layered and grouped according to the excavation steps to form a step-by-step excavation simulation model. The measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and the excavation calculation is carried out step by step. In each calculation step, when the stress of the joint unit reaches the yield condition, stress pull-back correction is performed based on the vertical partition of the plastic potential function, the inelastic deformation and joint mechanical parameters are iteratively updated, and the surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Based on the distribution of damage types, damage degree, and displacement vectors, the V-shaped failure zone of the top arch, the extrusion zone of the sidewalls, and the fault zone of the layered rock blocks are identified, and the large deformation analysis results of the deep layered fractured rock are obtained.

2. The method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to claim 1, characterized in that, The process involves acquiring geological survey data and 3D laser scanning point cloud data of deep layered fractured rock tunnels, and extracting geometric features of rock mass structural surfaces from the 3D laser scanning point cloud data to obtain a structural surface geometric feature dataset. Specifically: Geological survey data of deep layered fractured rock tunnels are obtained, stratigraphic lithology information, joint statistics and geological structure distribution information of the target area are extracted, and three-dimensional laser scanning equipment is used to perform full-area scanning of the tunnel face and exposed rock surface to obtain point cloud data of the rock surface. The point cloud data is subjected to noise reduction filtering and coordinate registration to construct a three-dimensional digital model of the rock mass surface. The three-dimensional digital model of the rock mass surface is then fitted with a local surface to calculate the normal vector of each data point in the point cloud. Based on the abrupt change feature of normal vector change, the structural surface traces of the rock mass are identified. The continuous structural surface boundary is traced along the trace direction. The closed structural surface boundary is divided into regions to obtain planar fragments of individual structural surfaces. Spatial plane fitting is performed on each planar fragment to calculate the spatial dip and dip angle of the planar fragment. The spacing between adjacent structural surface traces in the same plane is statistically analyzed. Multiple groups of planar fragments are grouped according to clustering thresholds of dip and tilt angle. The average spacing of structural surfaces and the number of structural surfaces per unit area in each group are calculated to obtain the density parameters of each group of structural surfaces. The spatial geometric parameters of structural surfaces in all groups are summarized to construct a structural surface geometric feature dataset containing structural surface dip, tilt angle, spacing, density, and distribution range.

3. The method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to claim 1, characterized in that, The discrete fracture network algorithm is used to reconstruct the geometric feature dataset of the structural surface to obtain a discrete fracture network model. Based on this model, equivalent mechanical parameters of rock masses of different sizes are calculated. When these parameters tend to a preset stable value as the rock mass size increases, the size of the rock mass characterizing unit is determined. Specifically: Using the dip, tilt angle, spacing and density parameters in the geometric feature dataset of the structural surface as input variables, the Monte Carlo random sampling method is used to generate multiple sets of independent random structural surfaces in the three-dimensional spatial domain. The intersection operation between each set of random structural surfaces and the boundary of the spatial domain is performed to extract the trace distribution features of the structural surface on the boundary. Based on the trace distribution characteristics, a three-dimensional discrete fracture network model consistent with the orientation distribution of the structural surface in the target region is constructed. Based on the three-dimensional discrete fracture network model, multiple cubic rock mass models are set in order of increasing size. The discrete fracture network in each cubic rock mass model is spatially cut to obtain a set of fragmented rock blocks of the corresponding size. Apply uniaxial compression boundary conditions to each set of fractured rock blocks, calculate the equivalent compressive strength of the set of fractured rock blocks at each size, and plot the convergence curve of the equivalent compressive strength as the size of the cubic rock mass model changes. When the change in equivalent compressive strength during the increase of the size of the cubic rock mass model is lower than the set convergence threshold, the size is determined as the size of the rock mass characterization unit.

4. The method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to claim 1, characterized in that, The method involves constructing a synthetic rock mass numerical model with elastic boundaries using the unit size of the rock mass characterizing element, coupling it with a three-dimensional bonded block damage constitutive model, defining the linear degradation law of joint normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle with inelastic deformation, and grouping the surrounding rock units into layers according to the excavation sequence to form a step-by-step excavation simulation model. Specifically: The analysis range of the discrete fracture network is defined based on the size of the rock mass characterization unit. All structural surfaces in the discrete fracture network are retained within the analysis range. The structural surfaces are used as cutting boundaries to spatially divide the rock blocks, resulting in a discretely distributed set of fractured rock blocks. Elastic boundary blocks with a size larger than the analysis range are generated at the critical periphery of the analysis range. The mechanical parameters of the elastic boundary blocks are consistent with the mechanical parameters of the intact rock blocks to eliminate the influence of boundary effects on the calculation results. The tunnel space is cut out inside the fragmented rock block set according to the design cross-sectional shape and size of the tunnel. The surrounding rock units within a predetermined range around the tunnel are divided into multiple excavation groups according to the actual excavation advance and construction sequence. Each excavation group corresponds to an independent excavation stage. The constitutive model of the intact rock block in the set of fractured rock blocks is set as an elastic constitutive model, and the constitutive model of the structural surface unit is set as a three-dimensional bonded block damage constitutive model. The three-dimensional bonded block damage constitutive model decomposes the displacement vector of the structural surface into elastic deformation part, plastic deformation part and fracture deformation part, and merges the plastic deformation part and fracture deformation part into inelastic deformation amount under tension, tension-shear and compression-shear stress states. Once the structural surface reaches the yield condition, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle are defined as linear functions of the inelastic deformation. Under compression, the structural surface exhibits only linear elastic behavior and does not enter the yield stage. Under tension and tension-shear conditions, the normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle all decrease linearly with the increase of inelastic deformation. Under compression-shear conditions, only tensile strength, cohesion, and dilatation angle decrease linearly with the increase of inelastic deformation. The normal stiffness and tangential stiffness remain unchanged from their initial values. After constitutive coupling and element grouping of the model, a step-by-step excavation simulation model containing elastic boundaries, fractured rock blocks, structural surfaces, and layered excavation elements is obtained.

5. The method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to claim 1, characterized in that, The measured in-situ stress boundary conditions are applied to the step-by-step excavation simulation model, and excavation calculations are carried out step by step. In each calculation step, when the stress of the joint element reaches the yield condition, stress pullback correction is performed based on the vertical partition of the plastic potential function, and the inelastic deformation and joint mechanical parameters are iteratively updated. The surrounding rock displacement field and structural surface damage type distribution are output simultaneously. Specifically: An initial geostress field measured on-site is applied to the step-by-step excavation simulation model to constrain the normal displacement of the model boundary. The empty constitutive structure of the corresponding surrounding rock unit is activated in the order of the excavation group to simulate the excavation unloading effect. In each calculation step, the maximum unbalanced force of the model tends to stabilize as the convergence criterion. Calculate the elastic test stress of the joint element, substitute the elastic test stress into the yield surface function for judgment, when the elastic test stress meets the yield condition, construct the plastic potential function under the current calculation step, the plastic potential function is the linear relationship between tangential stress and normal stress, draw a perpendicular line to the plastic potential function through the origin, and use the perpendicular line as the partition boundary to divide the stress region above the yield surface into the lower left region and the upper right region. For the elastic test stress of the joint element located in the lower left region, calculate the slope of the line connecting the elastic test stress and the origin, project the elastic test stress onto the yield surface along the direction of the line to obtain the corrected joint element stress, and extract the components of the direction of the line in the normal and tangential directions as inelastic strain components. For the elastic test stress of the joint element located in the upper right region, calculate the normal vector of the plastic potential function at the origin, project the elastic test stress onto the yield surface along the direction of the normal vector to obtain the corrected joint element stress, and calculate the inelastic strain components based on the components of the normal vector in the normal and tangential directions. Divide the inelastic strain components by the initial normal stiffness and the initial tangential stiffness to obtain the inelastic strain coefficients in the normal and tangential directions, respectively. Combine the corrected joint element stress to calculate the inelastic multiplier. Multiply the inelastic multiplier by the inelastic strain coefficients in the normal and tangential directions respectively and sum them to obtain the inelastic deformation increment of the current calculation step. Add the inelastic deformation amount of the previous calculation step to obtain the total inelastic deformation amount of the current calculation step. The normal stiffness, tangential stiffness, tensile strength, cohesion, and dilatation angle of the joint are updated according to the total inelastic deformation in accordance with the linear deterioration law. The damage degree of the bonded block is calculated based on the ratio of the updated normal stiffness to the initial normal stiffness. The displacement values ​​of each node of the surrounding rock are obtained by solving the overall equilibrium equation, and the displacement field of the surrounding rock is constructed. The damage type of the structural surface is determined according to the stress state and inelastic deformation component of the modified joint element. The damage type includes tensile damage, compression-shear damage and tension-shear damage. The surrounding rock displacement field, the distribution of structural surface damage types and the corresponding damage degree data are output simultaneously.

6. The method for analyzing large deformations of deep layered fractured rocks based on multi-source fusion according to claim 1, characterized in that, Based on the damage type distribution, damage degree, and displacement vector, the V-shaped failure zone of the crown arch, the extrusion zone of the sidewalls, and the faulting zone of layered rock blocks were identified, resulting in the analysis results of large deformation of deep layered fractured rock. Specifically: The displacement vector of each node in the surrounding rock displacement field output by each calculation step of the step-by-step excavation simulation model is retrieved. The damage type and corresponding damage degree of the joint unit at the location of each node are associated. The arch area and the sidewall area are divided with the center line of the roadway as the reference. The damage type of each joint unit in the arch area is traversed. When the damage type of the joint unit is tensile damage and the damage degree reaches the maximum value, and the displacement vector of the node corresponding to the joint unit points to the free face of the roadway, the area where the joint unit is located is marked as the arch V-shaped failure zone. The damage type of each joint unit in the sidewall region is traversed. When the damage type of the joint unit is compressive-shear damage and the damage degree reaches the maximum value, and the displacement vector of the corresponding node of the joint unit points to the free surface inside the roadway, the area where the joint unit is located is marked as the sidewall extrusion zone. Paired tracking of joint units on both sides of the layered structural plane. When both paired joint units exhibit tensile-shear damage and the damage degree continues to increase with the time step, and the displacement vectors of the corresponding nodes of the two joint units have relative shear displacement along the direction of the structural plane, the area where the structural plane is located is marked as the layered rock block displacement zone. By summarizing the spatial distribution range, corresponding damage extreme values, and displacement vector amplitudes of all marked V-shaped failure zones of the top arch, extrusion zones of the sidewalls, and fault zones of layered rock blocks, a mapping relationship between large deformation modes and spatial locations is constructed, and the large deformation analysis results of deep layered fractured rocks are obtained.

7. A system for analyzing large deformations of deep layered fractured rocks based on multi-source fusion, characterized in that, The deep layered fractured rock large deformation analysis system based on multi-source fusion includes a memory and a processor. The memory includes a program for a method of analyzing large deformation of deep layered fractured rocks based on multi-source fusion. When the program for analyzing large deformation of deep layered fractured rocks based on multi-source fusion is executed by the processor, it implements the steps of the method of analyzing large deformation of deep layered fractured rocks based on multi-source fusion as described in any one of claims 1 to 6.

Images