Spatial registration method for tunnel multi-source heterogeneous data

Abstract

The invention discloses a spatial registration method for tunnel multi-source heterogeneous data, and the method comprises the steps: constructing a lightweight digital twinborn model based on a design parameter library, a construction log library and a real-time monitoring library; mapping the multi-source sensor data in the real-time monitoring library to a digital twin model space node through a space-time synchronization algorithm; probability distribution model-based abnormal value elimination is performed on the geological exploration data, the while-drilling parameters and the geological sketch, and a discrete smooth interpolation algorithm is adopted to generate a surrounding rock mechanical parameter three-dimensional distribution field; the surrounding rock mechanical parameter three-dimensional distribution field serves as input, and front surrounding rock classification probability distribution is dynamically predicted through a machine learning model optimized through transfer learning; based on the three-dimensional point cloud data of the tunnel face in the digital twin model, extracting geometric parameters of the structural plane and constructing a fracture network model constrained by a normal vector; according to the method, the surrounding rock risk location is accurately positioned based on double-threshold triggering, and multi-scale model seamless fusion is realized through gradient constraint topology fusion.

Description

Technical Field

[0001] This application relates to the technical field of tunnel engineering, and particularly to a spatial registration method for multi-source heterogeneous data in tunnels. Background Art

[0002] With the rapid development of underground engineering construction, the geological conditions faced by tunnel construction are becoming increasingly complex, and the types of data involved in the construction process are numerous and heterogeneous in source, including geological exploration data (such as borehole cores, TSP seismic waves, ground penetrating radar images), construction monitoring data (such as mechanical sensor parameters, surrounding rock deformation monitoring values), equipment operation data (such as the spraying volume of wet shotcreting robots, positioning information of transport vehicles), etc. These data have significant differences in format, accuracy, and spatio-temporal scale, resulting in the following prominent problems for traditional technologies:

[0003] 1. The prediction of surrounding rock stability (such as classification probability) and the analysis of fracture networks are independent of each other, and only separated statistical indicators can be output (such as "80% probability of grade III surrounding rock" or "fracture density of 0.5 fractures per square meter"), and the geological and mechanical coupling risks in specific locations cannot be identified (for example: "the probability of instability of grade III surrounding rock at the crown of pile number K3+215 increases sharply due to a fracture connectivity rate > 60%").

[0004] 2. The global geological model (such as the surrounding rock mechanical parameter field) and the local engineering model (such as the support structure in the risk area) have geometric misalignment and attribute jumps (such as discontinuous boundary stress fields in sub-models) during fusion due to heterogeneous data sources and different scales, resulting in decision-making distortion. Summary of the Invention

[0005] To solve the above problems, an embodiment of the present invention provides a spatial registration method for multi-source heterogeneous data in tunnels, and the method includes:

[0006] Construct a lightweight digital twin model based on the design parameter library, construction log library, and real-time monitoring library;

[0007] Map the multi-source sensor data in the real-time monitoring library to the spatial nodes of the digital twin model through a spatio-temporal synchronization algorithm;

[0008] Perform outlier rejection on geological exploration data, drilling parameters, and geological sketches based on a probability distribution model, and use a discrete smooth interpolation algorithm to generate a three-dimensional distribution field of surrounding rock mechanical parameters;

[0009] Using the three-dimensional distribution field of surrounding rock mechanical parameters as input, dynamically predict the probability distribution of the surrounding rock grade ahead through a machine learning model optimized by transfer learning;

[0010] Extract the geometric parameters of the structural plane and construct a fracture network model with normal vector constraints based on the three-dimensional point cloud data of the tunnel face in the digital twin model;

[0011] When the probabilities of rock masses of grades I to III in the probability distribution of rock mass classification exceed the set threshold, and the fracture connectivity rate output by the fracture network model exceeds the critical value, the support parameters in the associated design parameter library, the geological event records in the construction log library, and the convergence deformation data are associated to construct a local digital twin model;

[0012] Perform multi-scale fusion of the local digital twin model and the three-dimensional distribution field of rock mass mechanical parameters based on topological similarity to eliminate geometric and property discontinuity regions.

[0013] Furthermore, for outlier rejection, a Gaussian mixture model is used to verify the consistency of data distribution, and outlier data deviating from the preset probability interval is rejected.

[0014] Furthermore, the machine learning model is a bidirectional LSTM network embedded with a multi-head attention mechanism, and the input parameters include the dynamic coupling characteristics of the rock mass integrity coefficient and the uniaxial compressive strength of the rock.

[0015] Furthermore, the normal vector constraint specifically includes: identifying the dominant fracture group based on normal vector space clustering, and the spatial autocorrelation characteristics of the fracture density distribution are fused in the clustering process.

[0016] Furthermore, for multi-scale fusion, a gradient constraint deformation algorithm is adopted, and the boundary deformation of the local digital twin model is adaptively adjusted according to the spatial gradient of the mechanical parameters of the surrounding rock in the neighborhood.

[0017] Furthermore, after identifying the dominant fracture group, a conditional generative adversarial network is used to construct the spatial extension probability field of the structural plane to predict the fracture development trend in the unexcavated area ahead.

[0018] Furthermore, when constructing the local digital twin model, the parametric BIM component library is called, and the component selection logic is based on the mapping relationship between geological anomaly types and the stability grade of the surrounding rock.

[0019] Furthermore, the spatio-temporal synchronization algorithm adopts multi-sensor fusion positioning technology to achieve sub-millimeter-level spatial positioning and millisecond-level time synchronization.

[0020] The technical effects and advantages of a spatial registration method for multi-source heterogeneous data of tunnels provided by the present invention:

[0021] The present invention accurately locates the risk areas of surrounding rock based on double-threshold triggering, and realizes seamless fusion of multi-scale models through gradient-constrained topological fusion; the present invention predicts the probability distribution of the surrounding rock grading ahead through a bidirectional LSTM network optimized by transfer learning, and uses a fracture network model constrained by normal vectors to calculate the connectivity rate in real time, thereby dynamically coupling the surrounding rock grading and the fracture network; when the probability of grade I-III surrounding rock > threshold and the fracture connectivity rate exceeds the critical value, the support parameters and geological event records are automatically associated to construct a locatable local sub-model, thereby accurately locating the risk area; the parametric BIM component library is used to match the geological anomaly type and the support structure to ensure smooth transition of mechanical properties at the fusion boundary; the mechanical field of the surrounding rock generated based on the discrete smooth interpolation algorithm constrains the topological similarity in the process of sub-model fusion and eliminates geometric distortion. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 is a flowchart of the spatial registration method for multi-source heterogeneous data of a tunnel in an embodiment of the present application;

[0023] Figure 2 is a schematic diagram of the real-time monitoring library of multi-source heterogeneous data of a tunnel in an embodiment of the present application;

[0024] Figure 3 is a schematic diagram of the digital twin model of the lithology of the surrounding rock in an embodiment of the present application;

[0025] Figure 4 is a schematic diagram of the coupling surrounding rock grading decision of the uniaxial compressive strength of rock and the rock mass integrity coefficient in an embodiment of the present application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0026] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0027] Please refer to Figure 1 as shown, the embodiment of the present invention provides a spatial registration method for multi-source heterogeneous data of a tunnel, and the method includes:

[0028] S1: Based on the design parameter library, construction log library and real-time monitoring library, construct a lightweight digital twin model;

[0029] S2: Map the multi-source sensor data in the real-time monitoring library to the spatial nodes of the digital twin model through a spatio-temporal synchronization algorithm;

[0030] S3: Remove outliers from geological exploration data, drilling parameters, and geological sketches based on a probability distribution model, and generate a three-dimensional distribution field of surrounding rock mechanical parameters using a discrete smooth interpolation algorithm;

[0031] S4: Using the three-dimensional distribution field of surrounding rock mechanical parameters as input, dynamically predict the probability distribution of the surrounding rock classification ahead through a machine learning model optimized by transfer learning;

[0032] S5: Based on the three-dimensional point cloud data of the tunnel face in the digital twin model, extract the geometric parameters of the structural plane and construct a fracture network model with normal vector constraints;

[0033] S6: When the probabilities of surrounding rock grades I to III in the probability distribution of surrounding rock classification exceed the set threshold, and the fracture connectivity rate output by the fracture network model exceeds the critical value, associate the support parameters in the design parameter library, the geological event records in the construction log library, and the convergence deformation data to construct a local digital twin sub-model;

[0034] S7: Perform multi-scale fusion based on topological similarity between the local digital twin sub-model and the three-dimensional distribution field of surrounding rock mechanical parameters to eliminate geometric and attribute discontinuity regions.

[0035] It should be specifically noted that in the present invention, the digital twin model has multiple variants, such as the digital twin model of surrounding rock moisture content, the digital twin model of surrounding rock lithology, and the digital twin model of surrounding rock integrity, etc.

[0036] In step S3, it is necessary to preprocess three types of heterogeneous data sources from geological exploration (such as core drilling tests), drilling monitoring (such as drilling speed, torque), and manual geological sketches. The core objective is to remove outliers to ensure the reliability of the subsequent generated three-dimensional distribution field of surrounding rock mechanical parameters (such as elastic modulus, compressive strength, etc.). These data have significant differences in acquisition methods, accuracy, density, and inherent errors, which will inevitably result in the existence of abnormal data points in the original dataset that deviate significantly from the main distribution and may be caused by measurement errors, local extreme geological bodies, or recording errors.

[0037] The specific method for removing outliers is as follows:

[0038] Align the discrete point mechanical parameters obtained from geological exploration (such as the uniaxial compressive strength at a depth of 10 m in a certain borehole), the drilling parameters continuously recorded during drilling (such as drilling speed, torque, usually sampled by time or depth), and the qualitative or semi-quantitative rock mass information described in geological sketches (such as "extremely developed joints", which need to be converted into corresponding quantitative index ranges) in space and time and unify the format to form a preliminarily integrated dataset; each data point contains its spatial position coordinates (X, Y, Z) and one or more corresponding mechanical parameter observation values;

[0039] Construct a Gaussian Mixture Model (GMM). For a single mechanical parameter of the surrounding rock to be processed (for example, in this processing, the focus is on the elastic modulus), extract all the observed values of this parameter from the preliminarily fused dataset;

[0040] Assume that within a relatively homogeneous geological unit, the observed values of this parameter may be generated by multiple different but inherent geological states (such as different lithologies, different weathering degree regions), and each state corresponds to a Gaussian distribution (normal distribution);

[0041] Use the Gaussian Mixture Model (GMM) to perform probability modeling on the set of observed values of the extracted elastic modulus. The GMM model regards the overall data distribution as the weighted sum of K (K needs to be determined according to the actual data characteristics and geological understanding) Gaussian distribution components. The GMM model automatically learns through algorithms such as Expectation-Maximization (EM), that is, the mean (representing the central tendency of this subclass of data), variance (representing the dispersion degree of this subclass of data), and the weight of this component in the entire mixture model (representing the proportion of this subclass of data) of each Gaussian component.

[0042] The trained GMM model provides the probability density distribution of the data space. For each observed value V of the elastic modulus in the dataset i , the model can calculate the probability density value P(V i ) that it belongs to this mixture model, or more commonly, calculate its log-likelihood value. The log-likelihood value reflects the possibility that V i conforms to the "normal" data distribution described by the GMM model;

[0043] Set a preset probability (or likelihood) interval as the judgment criterion. For example, a threshold based on the log-likelihood value distribution of all data points can be set (such as 3 times the standard deviation lower than the average log-likelihood value of all data points), or a very low probability density threshold can be set.

[0044] Traverse all the observed values of the elastic modulus. If the probability density P(V i ) or the log-likelihood value of a data point is significantly lower than the preset log-likelihood threshold P threshold , that is, P(V i ) < P threshold or its log-likelihood value is far lower than the normal range, then it is considered that this data point "deviates from the preset probability interval", is inconsistent with the main data distribution described by the GMM model, and is very likely to be an outlier. This data point will be marked and removed from the current processed elastic modulus dataset.

[0045] Exemplary:

[0046] Suppose in a specific section of a tunnel, based on geological understanding, the surrounding rock is mainly medium to slightly weathered granite, and its elastic modulus is expected to be mainly distributed between 10 GPa and 30 GPa; through GMM modeling, GMM modeling may learn 1 - 2 main distribution components (for example, one concentrated around 15 GPa and one around 25 GPa), as well as their respective variances and weights;

[0047] If the elastic modulus value recorded at a depth of 20 m in a certain borehole is 0.5 GPa (far lower than the expected lower limit), or the instantaneous elastic modulus inferred from a certain measurement - while - drilling point is as high as 100 GPa (far higher than the expected upper limit), the probability densities calculated by GMM for these will be extremely low (for example, far lower than the set P threshold ), and they will be identified as outliers and removed.

[0048] Similarly, if a data point falls within the main distribution range (such as 18 GPa), but the values of all adjacent points around its occurrence are within the range of 25 - 30 GPa (forming a local "cluster"), and this point itself is significantly lower than the sub - distribution formed by its adjacent points in the local probability density calculation of GMM, it may also be identified as an anomaly due to its "inconsistency" in this local context (if the model can capture this local pattern).

[0049] The above process usually needs to be carried out separately for each key surrounding - rock mechanical parameter (such as compressive strength, cohesion, internal friction angle, etc.). After removing outliers, the resulting dataset is cleaned and better reflects the distribution of the mechanical properties of the real geological body; these cleaned data will be used as the input for the subsequent Discrete Smooth Interpolation (DSI) algorithm to generate a continuous, smooth and physically - reasonable three - dimensional distribution field of the surrounding - rock mechanical parameters.

[0050] In step S4, it is necessary to dynamically predict the probability distribution of the surrounding - rock classification in the unexcavated area ahead of the tunnel using the three - dimensional distribution field of the surrounding - rock mechanical parameters (including spatially continuous attributes such as elastic modulus and compressive strength) (such as the probabilities of grades I - V in the "Standard for Engineering Rock Mass Classification"), and the machine - learning model architecture and input - feature construction methods include:

[0051] The architecture design of the machine - learning model includes:

[0052] The core uses a Bidirectional Long Short - Term Memory Network (Bi - LSTM) as the basic framework. Bi - LSTM can simultaneously learn the forward dependence (from historical points to the current point) and backward dependence (from future potential points to the current point) of the geomechanical parameters in the tunnel - driving direction (i.e., the spatial - sequence direction, such as along the tunnel axis from the excavated section to the unexcavated section), and effectively capture the long - distance correlation of the surrounding - rock parameters in the spatial sequence.

[0053] Above the Bi-LSTM network layer, a multi-head attention mechanism is embedded. The multi-head attention mechanism allows the model to concurrently focus on information at different positions, different feature dimensions, or different geologically significant segments in the input sequence, specifically including:

[0054] Use the hidden state sequence output by Bi-LSTM as the input to the multi-head attention layer;

[0055] The multi-head attention layer maps the input features to multiple subspaces (referred to as "heads"), and independently calculates the association weights (i.e., attention scores) between all positions in the sequence within each subspace.

[0056] Each "head" can focus on learning a specific dependency pattern (e.g., one head focuses on the mutation points of rock mass integrity, one head focuses on the uniaxial compressive strength of rock, and the last head focuses on the active groundwater area).

[0057] Finally, the attention results of all "heads" are concatenated and linearly transformed to form a weighted feature representation that integrates global context information, which significantly enhances the model's ability to identify key risk areas under complex geological conditions.

[0058] The main input of the machine learning model is the slice sequence data of the three-dimensional distribution field of the surrounding rock mechanical parameters generated in step S3 along the tunnel axis direction (driving direction); each slice represents the spatial distribution (gridded data) of the mechanical parameters at a certain distance in front of the tunnel face (e.g., one slice per meter or every 5 meters).

[0059] Introduce the dynamic coupling characteristics of the rock mass integrity coefficient (Kv) and the uniaxial compressive strength of rock (Rc) as additional inputs:

[0060] The rock mass integrity coefficient (Kv) is derived from the fracture network model constructed by extracting the structural plane parameters based on the three-dimensional point cloud data of the tunnel face. Kv quantitatively characterizes the degree of fragmentation of the rock mass cut by joints and fractures (the value range is 0 - 1, and the lower the value, the more fragmented); this parameter changes dynamically in space and is updated in real time as new point cloud data is obtained and the fracture model is updated in front of the driving face.

[0061] The uniaxial compressive strength of rock (Rc) is derived from the real-time monitoring data mapped to the nodes of the digital twin model such as Figure 3 through the spatio-temporal synchronization algorithm, which mainly reflects the ability of the rock itself to resist axial pressure failure, and its strength is affected by many factors (such as Kv and groundwater conditions);

[0062] As Figure 4 shown, the machine learning model does not simply take Kv and Rc as independent features for input, but constructs a coupling feature that can reflect the interaction between the two based on the fuzzy function and the three-dimensional point cloud data of the tunnel face;

[0063] By leveraging the sequence modeling ability of Bi-LSTM itself, the network learns the dynamic association patterns between the components of Kv and Rc in the temporal and spatial sequences (whether the continuous decrease of the forward Kv value (the rock mass is becoming more fragmented) is accompanied by a synchronous decrease in the Rc value (the strength of the rock itself is also deteriorating), or whether Rc remains relatively stable).

[0064] The multi-head attention mechanism further helps the model focus on the most predictive coupling relationships between Kv and Rc at specific positions (for example, it can help the model identify the abnormal points where Rc and Kv are "mismatched", such as areas with high Rc (hard rock) but extremely low Kv (extremely fragmented) (high-risk areas for rock bursts), or areas with medium Rc but extremely high Kv (possibly intact massive layered rock masses)).

[0065] The training method of the machine learning model includes;

[0066] For the current face position, the forward prediction area is divided into continuous prediction units (slices); the input of the machine learning model is a spatial sequence along the axis direction, from the excavated section (with known true surrounding rock stability grades) to a certain distance behind the current face (providing context), and then to multiple units to be predicted in front of the face (each unit contains statistical features of the mechanical parameter field, coupling features of Kv and Rc, etc.) of this slice.

[0067] The input of the machine learning model includes: the machine learning model outputs a probability distribution vector (such as [P(Level Ⅰ), P(Level Ⅱ), P(Level Ⅲ), P(Level Ⅳ), P(Level Ⅴ)]) for each forward prediction unit, indicating the probability that the unit belongs to different surrounding rock stability grades.

[0068] Use the historical data of other tunnel projects with similar geological conditions in the integrated design parameter library in step S1 (including complete mechanical parameter fields, Kv, Rc records, and finally verified surrounding rock stability grades) to pre-train the machine learning model; when applied to this tunnel, use the real-time data accumulated in the current tunnel (new known points are continuously generated as tunneling progresses) to fine-tune the pre-trained model; this effectively solves the problem of insufficient initial data for new tunnels and improves the generalization ability and prediction accuracy of the model.

[0069] As the tunnel continues to be excavated, new face point cloud data (updating Kv), new monitoring data (updating Rc), new mechanical parameters (the distribution field is updated after the drilling data is processed by S3), and the information of the actually exposed surrounding rock stability grades are continuously added. The input sequence of the machine learning model is updated accordingly, and online fine-tuning can be performed periodically to achieve dynamic and progressive prediction optimization.

[0070] Exemplary:

[0071] Assume the tunnel is driven to station K10+500; the input sequence of the machine learning model includes:

[0072] From K10+450 to K10+499 (excavated section, with known true surrounding rock stability grades, used to provide context and for training);

[0073] K10+500 (current tunnel face, with partial information known);

[0074] From K10+501 to K10+550 (section to be predicted);

[0075] For the element K10+520 in the predicted section:

[0076] The input includes the mechanical parameters of the slice at this location (such as average elastic modulus of 25 GPa and standard deviation of 3 GPa).

[0077] The input includes the rock mass integrity coefficient Kv = 0.65 (medium integrity) at this location (predicted or interpolated based on the forward point cloud).

[0078] The input includes the real-time monitored interpolation / prediction of the groundwater state at this location: water inflow = 15 L / min, pore water pressure = 200 kPa.

[0079] The dynamically coupled features constructed by the model may include Kv * water inflow = 0.65 * 15 = 9.75, and pore water pressure / Kv = 200 / 0.65 ≈ 307.7.

[0080] After being processed by the Bi-LSTM and multi-head attention mechanism, the machine learning model may output the surrounding rock classification probabilities for this element as: P(Grade Ⅲ) = 60%, P(Grade Ⅳ) = 35%, P(Grade Ⅴ) = 5%. This indicates that there is a certain risk at this location (the combined probability of Grade Ⅳ and above is 40%), and close attention is required; this prediction result will be used for the threshold judgment in step S6.

[0081] In step S5, it is necessary to extract the geometric parameters (such as position, attitude, trace length, aperture) of the structural planes (joints, fractures, bedding planes, etc.) based on the 3D point cloud data of the tunnel face, and finally construct a fracture network model constrained by the normal vector; the core method of normal vector constraint, that is, to identify the dominant fracture groups by fusing the spatial autocorrelation characteristics of the fracture density distribution through normal vector space clustering:

[0082] Data basis and preprocessing include:

[0083] The input data is the 3D point cloud of the current tunnel face obtained by a laser scanner or photogrammetry; this point cloud densely represents the geometric morphology of the tunnel face surface and contains a large amount of discontinuous surface information formed by the structural planes (fractures).

[0084] Identify each independent structural patch from the point cloud using point cloud segmentation and plane fitting algorithms (such as RANSAC, region growing method); operations for each identified structural patch include:

[0085] Calculate the spatial coordinates (X, Y, Z) of its center point;

[0086] Calculate its normal vector (Nx, Ny, Nz) by fitting a plane. The normal vector is a unit vector perpendicular to the plane of the structural surface, uniquely defining the spatial orientation of the structural surface (dip direction and dip angle can be obtained by converting the normal vector);

[0087] Extract or estimate other geometric parameters, such as trace length (the visible length of the structural surface on the outcrop), aperture (the degree of crack opening), etc.

[0088] Traditional clustering (such as K-means or DBSCAN based solely on the direction of the normal vector) may ignore the clustering characteristics of the spatial distribution of structural surfaces (for example, a certain set of cracks may be particularly developed in a certain section of the tunnel, forming a dense zone); the identification of dominant crack groups by normal vector spatial clustering includes: classifying structural surfaces with similar spatial orientations (i.e., similar normal vector directions) into the same dominant crack group. The dominant crack group represents the main structural surface directions that control the engineering mechanical behavior of the rock mass (such as strength, permeability, failure mode);

[0089] Clustering methods that incorporate spatial autocorrelation include:

[0090] For each structural patch i, construct a composite feature vector that includes its spatial position information and normal vector information. For example, Feature i = [X i , Y i , Z i , Nx i , Ny i , Nz i , where (X i , Y i , Z i ) is the coordinate of its center point, and Nx i , Ny i , Nz i is its unit normal vector.

[0091] Methods for calculating spatial autocorrelation weights include:

[0092] Divide the study area into regular grid cells in three-dimensional space (or in a two-dimensional section projected onto the tunnel axis direction);

[0093] Statistical crack density within each grid cell (such as the number of cracks or total trace length per unit area or unit volume);

[0094] Calculate the spatial autocorrelation of fracture density, which measures the similarity of fracture density in adjacent regions in space. For example, if the fracture density of unit A is high and the fracture density of its adjacent unit B is also high, then A and B have positive spatial autocorrelation (aggregation) in terms of fracture density; if A has a high density while B has a low density, it shows negative spatial autocorrelation (dispersion).

[0095] Based on the results of spatial autocorrelation weight analysis (such as confirming the existence of autocorrelation by calculating the global or local Moran's I index), define a spatial density similarity weight W ij (spatial) for each pair of structural plane patches (i, j); this weight reflects the strength of the spatial correlation of the fracture density at the locations of structural planes i and j; if i and j are located in a region with high positive autocorrelation of fracture density (i.e., they are both in a dense zone or both in a sparse zone), then the value of W ij (spatial) is large; if they are located in a dense zone and a sparse zone respectively (negative correlation), then the value of W ij (spatial) is small or negative.

[0096] Adopt a clustering algorithm that can incorporate custom distance or similarity metrics (such as spectral clustering, DBSCAN with custom distance), and the methods for defining the comprehensive distance and similarity between structural planes i and j include:

[0097] Calculate the cosine of the angle between the normal vectors of i and j, Cosθ = N i ·N j (dot product, since they are unit vectors). The closer Cosθ is to 1, the more similar the directions are; at the same time, the obtained W ij (spatial);

[0098] Combine the cosine of the angle with the spatial density similarity weight. For example, Similarity ij =α*|Cosθ| + β*W ij (spatial); where α and β are weight coefficients that adjust the importance of the two (which can be determined by geological experience or experiments), and |Cosθ| ensures that only the magnitude of the direction angle (0° - 90°) is concerned, without distinguishing

[0099] between positive and negative normal vectors (representing the same plane).

[0100] Use the above comprehensive similarity Similarity ij as the input and execute the clustering algorithm (such as setting a similarity threshold and the minimum number of cluster points), which will cluster the structural planes with similar normal vector directions and similar spatial distribution patterns (density autocorrelation) into clusters.

[0101] Finally, multiple dominant fracture groups are output, and each dominant fracture group contains a set of structural planes. These structural planes not only have similar spatial orientations but also often show a consistent pattern of aggregation or dispersion in space (guaranteed by the fused spatial autocorrelation characteristics).

[0102] The method for constructing a fracture network model with normal vector constraints includes:

[0103] For each identified dominant fracture group, calculate the average direction of the normal vectors of all the structural planes in this dominant fracture group (or use the representative normal vector of this group);

[0104] Statistically analyze the distribution of other geometric parameters of the structural planes in this group (such as the probability distribution models of trace length, aperture, and spacing);

[0105] Based on the above average direction and spatial position information, use the stochastic discrete fracture network (DFN) modeling technique, that is:

[0106] In the three-dimensional geological model space, according to the average normal vector direction (i.e., attitude) of each dominant fracture group, the statistical distribution of geometric parameters, and the density pattern of the actual distribution of the structural planes in this group in space (constraining the distribution range or intensity by the results of the spatial autocorrelation analysis in step S2), randomly generate a large number of virtual fractures that conform to the statistical laws of this dominant fracture group;

[0107] The attitude (normal vector direction) of the generated virtual fractures strictly follows the statistical characteristics (mean value, variance) of the dominant fracture group to which they belong, rather than being completely random, which ensures that the generated fracture network is highly consistent with the actual exposure situation in the orientation of the main structural planes;

[0108] Overlay the virtual fracture networks generated by all dominant fracture groups, that is, form the final fracture network model with normal vector constraints for engineering analysis; this fracture network model can more realistically reflect the spatial distribution and geometric characteristics of the key fracture groups that control the rock mass structure.

[0109] Exemplary:

[0110] Suppose 200 structural plane patches are identified in the point cloud of a certain tunnel face; traditionally, simply clustering based on normal vectors may result in 3 groups (for example: Group A: dip direction 120°, dip angle 70°; Group B: dip direction 30°, dip angle 50°; Group C: dip direction 210°, dip angle 80°).

[0111] Spatial autocorrelation analysis finds that the fracture density in the left area of the tunnel (from station K10 + 500 to K10 + 520) is abnormally high, and this high-density area has strong positive spatial autocorrelation (that is, the dense areas appear in patches).

[0112] In this area, the number of structural planes in both Group A and Group B has increased significantly and is distributed in a mixed manner.

[0113] After fusing spatial weight clustering:

[0114] Group A is subdivided into two subgroups: A1 (mainly distributed in the high-density area on the left, with an average trend of 125° and an inclination angle of 72°) and A2 (sparsely distributed in other areas, with an average trend of 115° and an inclination angle of 68°).

[0115] Group B is also subdivided into B1 (high-density area on the left, with an average trend of 35° and an inclination angle of 48°) and B2 (other areas, with an average trend of 25° and an inclination angle of 52°).

[0116] There is no obvious spatial aggregation in Group C, and it remains as one group.

[0117] When constructing the fracture network model:

[0118] In the high-density area on the left (K10+500-520), virtual fractures will be created according to the average normal vector direction, geometric parameter statistical values of A1 and B1, and a higher generation density.

[0119] In other areas, it is generated according to the statistical attributes of Group A2, B2, and C and a lower density.

[0120] The results include: The generated fracture network model can not only reflect the main fracture directions, but also accurately depict the concentrated development characteristics of multiple groups of dominant fractures in the key risk area (high-density area on the left); this will make the fracture connectivity rate calculated based on this model in step S6 (especially the connectivity within the high-density area) more accurate and reliable.

[0121] In step S7, it is necessary to seamlessly fuse the local digital twin sub-model (such as a higher-precision mechanical-seepage coupling model established for high-risk fault fracture zones, water-rich sacs, or large deformation areas) with the overall tunnel digital twin model (including the probability distribution of surrounding rock classification, the distribution field of mechanical parameters, the fracture network model, etc.). That is, the multi-scale fusion method includes using a gradient-constrained deformation algorithm to adaptively adjust the boundary deformation of the local digital twin sub-model according to the spatial gradient of the mechanical parameters of the adjacent surrounding rock.

[0122] A three-dimensional grid model (such as a tetrahedral or hexahedral grid) covering the entire tunnel engineering scope, where the grid nodes carry the three-dimensional distribution field of the mechanical parameters of the surrounding rock generated in step S3 (such as elastic modulus, compressive strength), the probability distribution of the surrounding rock classification predicted in step S4, the equivalent permeability and strength parameter fields derived from the normal vector-constrained fracture network model constructed in step S5, etc. This three-dimensional grid model has a relatively low resolution (such as a grid size of several meters to ten meters), and is used for overall stability analysis and macroscopic decision-making.

[0123] The local digital twin model includes: a high-resolution refined model (with a grid size up to the centimeter to decimeter level) established for the identified high-risk areas (such as the section from station number K10+520 to K10+540, where the predicted probability of grade IV is >60% and the fracture connectivity rate > the critical value); the local digital twin model adopts more complex constitutive relations (such as elastoplasticity, fluid-solid coupling) to accurately simulate the stress, strain, seepage response and potential failure mode of this area under excavation disturbance.

[0124] The mechanical properties (such as elastic modulus) of the boundary region of the local digital twin model are usually spatially gradual or have mutations in the global model. Traditional rigid splicing or simple interpolation will lead to discontinuous stress and deformation at the boundary, destroying the physical rationality of the solution.

[0125] The fusion requirements include: the calculation results (such as deformation, stress, plastic zone) of the local digital twin model need to be inversely mapped back to the global model to update the twin state; at the same time, the global model needs to provide constraints that conform to the actual geological conditions for the boundary of the local digital twin model.

[0126] The gradient constraint deformation algorithm includes:

[0127] Define an annular overlapping transition zone with a certain thickness between the local digital twin model and its surrounding global model. The annular overlapping transition zone is simultaneously covered by the tunnel's overall digital twin model grid and the local digital twin model grid;

[0128] In the annular overlapping transition zone, select the key surrounding rock mechanical parameters (such as elastic modulus (E)) on the grid nodes of the tunnel's overall digital twin model as the basis for gradient calculation;

[0129] Calculate the spatial gradient vector ∇E=(∂E / ∂x, ∂E / ∂y, ∂E / ∂z) of this parameter at each global grid node; this vector points to the direction in which the parameter value increases fastest, and its modulus |∇E| characterizes the severity of the parameter change (i.e., spatial variability);

[0130] According to the calculated spatial gradient information, construct a weight field W(x, y, z). This weight field is positively correlated with the modulus of the gradient |∇E|. For example, W(x, y, z)=|∇E| / max(|∇E|) or W(x, y, z)= sigmoid(|∇E|); thus, in the area where the mechanical parameters change violently (high-gradient area, such as fault interfaces, lithological boundaries), the weight is large; in the area where the parameters change gently (low-gradient area, such as inside homogeneous rock masses), the weight is small; this weight field is defined on the grid nodes of the tunnel's overall digital twin model in the annular overlapping transition zone.

[0131] During the solution process of the local digital twin model, the displacement degrees of freedom of its boundary nodes (located in the annular overlapping transition zone) are not completely fixed, but are affected by the constrained displacements provided by the global model at this position; when solving the local digital twin model, its boundary conditions are set as these displacement constraints adjusted by gradient weights.

[0132] The local digital twin model completes high-precision solution under the adjusted boundary conditions, and obtains high-resolution results (such as displacements, stresses, etc.) in the overlapping transition zone and inside.

[0133] The high-resolution results obtained by solving the local digital twin model (especially the results in the annular overlapping transition zone) are used to update the state variables (such as displacement field, stress field, plastic state flag) of the tunnel's overall digital twin model at the grid nodes in the corresponding area through conservation mapping or weighted average technology.

[0134] The updating method of the local digital twin model includes: using the solution of the updated tunnel's overall digital twin model to recalculate the gradient weights, then adjusting the boundary constraints of the local digital twin model, solving the local digital twin model, and backpropagating the results until the solution in the overlapping transition zone reaches a satisfactory consistency.

[0135] In the area where the spatial variation of the mechanical properties of the surrounding rock is gentle (low |∇E|, small W), the boundary deformation of the local digital twin model is mainly constrained by the global model to maintain overall coordination; the algorithm allows the boundary of the local digital twin model to produce greater adaptive deformation, enabling it to more accurately simulate the concentrated deformation or stress redistribution in local high-risk areas (such as weak interlayers), and avoiding false stress concentration or discontinuity.

[0136] Exemplarily:

[0137] Suppose there is a high-gradient area near the pile number K10+530: The global model shows that the elastic modulus E drops steeply from 30 GPa of hard granite to 0.5 GPa of fault gouge within a short distance (|∇E| is extremely large, W≈1).

[0138] A node I on the boundary of the local digital twin model (which finely simulates this fault zone).

[0139] The original solution of the global model gives the displacement of this point as 5 mm.

[0140] The average displacement of the high-resolution nodes adjacent to the inside of the local digital twin model is 20 mm (reflecting the large deformation trend of the fault gouge).

[0141] The local digital twin model is solved under this adjusted boundary displacement constraint. The actual calculated displacement of boundary node I may be close to 17 mm, significantly higher than the original 5 mm constraint; this allows the fault zone to produce the due large compressive deformation.

[0142] Meanwhile, in the low-gradient homogeneous rock mass area far from the fault (|∇E|≈0, W≈0), the boundary displacement of the local digital twin model basically remains unchanged (such as 8 mm), ensuring coordination with the overall model.

[0143] Finally, the large deformation result (17 mm) of the fault zone is updated to the nodes at this position in the global model through mapping, enabling the digital twin to truly reflect the actual state of the local high-risk area.

[0144] Furthermore, the existing DFN model can only characterize the fracture structure in the currently exposed area (the tunnel face and its vicinity), and cannot reliably predict the development trend of each dominant fracture group (such as the extension direction, density change, and enhanced connectivity) in the unexcavated area in front of the tunnel (especially in the high-risk section). Therefore, a generative model that can learn the spatial distribution law of geological structural planes needs to be constructed to output the spatial extension probability field of the structural planes. The spatial extension probability field of the structural planes quantitatively characterizes the probability that a structural plane of a specific dominant fracture group appears or the probability that its geometric parameters (such as density) reach a certain level at any position in the unexcavated area in front, providing a key basis for early warning (step S6).

[0145] The basic inputs of the CGAN include: based on the current tunnel face point cloud data, multiple dominant fracture groups (each group contains a set of structural planes with similar normal vector directions and consistent spatial distribution patterns); the statistical characteristics of the geometric parameters of each dominant fracture group (such as the mean / variance of the trace length, aperture distribution, spatial density); the normal vector constrained fracture network model (DFN) constructed by fusing spatial autocorrelation.

[0146] The architecture and training of the conditional generative adversarial network (CGAN) include:

[0147] The network structure includes a generator and a discriminator:

[0148] The generator includes: an input random noise vector (providing generation diversity) and a conditional information vector. The generator outputs a three-dimensional probability field grid, that is, a fake probability field. The value of each voxel in the grid represents the probability that the position belongs to the structural plane of a specific dominant fracture group or the probability of the fracture density level.

[0149] The discriminator includes: inputting real three-dimensional fracture distribution data (from the DFN model of a known area or high-precision detection) or the false probability field output by the generator; at the same time, the discriminator also receives the same conditional information vector; the goal of the discriminator is to distinguish whether the input data is "real" or "generated" and evaluate its matching degree with the conditional information; the conditional information vector is the key to the effective learning of CGAN, which encodes the geological background and known information that control the fracture extension.

[0150] Training includes:

[0151] Source of training samples: Utilize the complete geological data of other completed tunnel projects in the integrated design parameter library in step S1;

[0152] For example, a certain known section in the complete geological data;

[0153] Input includes (conditions): Take the characteristics of the dominant fracture set, geological structure background, etc. of the first half of this section (simulating "excavated and exposed") to construct a conditional vector.

[0154] Output includes (true label): Based on the actual exposure or fine detection data of the second half of this section (simulating "unexcavated"), construct a real three-dimensional fracture presence probability field or density field (for example, for grid cells, mark 1 if there is a fracture, otherwise mark 0; or calculate the total trace length of fractures within the cell as the density).

[0155] Adversarial training includes: Training CGAN with a large number of such samples; the generator learns to generate a realistic probability field that is sufficient to "deceive" the discriminator under a given conditional vector.

[0156] The discriminator learns to distinguish between the real probability field Real and the generated field Fake, and at the same time determines whether both are consistent with the conditional vector.

[0157] The training goal is to reach the Nash equilibrium, that is, the generator can generate a probability field that is almost indistinguishable from the real data distribution and conforms to the conditional vector of the geological conditions; at the same time, the discriminator cannot reliably distinguish between true and false.

[0158] Exemplary illustration:

[0159] Scenario: At the tunnel mileage K10+500, a set of steeply inclined joint sets (dominant set A) is exposed, with the average normal vector showing a dip of 210° and a dip angle of 75°; the advanced geological prediction shows that there is a small fault (strike N40°E) near K10+510 ahead.

[0160] Construction of the conditional vector (for set A):

[0161] Characteristics of set A: [Nx=-0.21, Ny=-0.58, Nz=0.79, L mean= 3.2 m, σL = 1.1 m, ρ current = 0.8 fractures / m²], where L mean is the average trace length, σL is the standard deviation of the trace length, and ρ current is the areal density of the current exposed area.

[0162] Geological structure: Fault location (K10 + 512), strike (N40°E), dip SE, dip angle 60°.

[0163] Current spatial pattern: High density in the exposed area and an extension trend along the 210° direction.

[0164] Encode to form the condition vector C A .

[0165] CGAN prediction:

[0166] Input C A and the noise z into the pre-trained generator G.

[0167] Output P A(X，Y，Z) : Predict the fracture presence probability field of Group A in the area from K10 + 501 to K10 + 550.

[0168] Result interpretation:

[0169] In the fault influence zone (K10 + 508 to K10 + 516), the P A value is generally > 0.85 (dark red), indicating that this group of fractures is very likely to be highly developed and have enhanced connectivity here (affected by fault disturbance).

[0170] Beyond K10 + 530, the P A value drops to 0.3 - 0.5 (light yellow), indicating that this group of fractures may gradually become sparse or pinch out.

[0171] Input the prediction result of this high-probability development area (K10 + 508 to K10 + 516) into step S6; combining the predicted low surrounding rock stability grade (from S4) and high water inrush risk (from S2) at this location, the system will trigger a high-level warning, prompting to strengthen the support and drainage measures in this section.

[0172] In step S6, when the surrounding rock classification probability distribution (from S4) and the normal vector constrained fracture network model (from S5) identify local risk areas (such as high-risk faults, water-rich pockets, large deformation sections), it is necessary to construct a refined mechanical-seepage local digital twin sub-model for special safety assessment; in this embodiment, by calling the parametric BIM component library and performing intelligent component selection according to the mapping relationship between geological anomaly types and surrounding rock stability grades, the efficient and standardized construction of the local digital twin sub-model is realized;

[0173] The BIM (Building Information Modeling) component library is a collection of digital models of predefined and parametric structural components for tunnel engineering; each component not only includes geometric shapes but also embeds key physical properties (such as material strength, permeability coefficient, support resistance) and behavioral rules (such as the constitutive relationship with surrounding rock).

[0174] The BIM component library is hierarchically organized according to functional types and applicable geological scenarios, including:

[0175] First-level classification (function): Initial support components (such as shotcrete layer, steel arch, bolt / cable system), advanced support components (such as pipe shed, advanced small duct, curtain grouting body), drainage components (such as drainage holes, blind drainage pipes, waterproof board), temporary reinforcement components (such as temporary invert, cross bracing).

[0176] Second-level classification (geological scenario): Under each type of functional component, it is further subdivided according to the dominant geological anomaly type and the range of surrounding rock stability grades applicable to its design;

[0177] Exemplary:

[0178] The shotcrete layer is further subdivided into: homogeneous low-risk type, fractured zone reinforcement type, high water inrush pressure type;

[0179] The bolt system is further subdivided into: short and dense type (soft surrounding rock), long and strong type (blocky surrounding rock), pressure dispersion type (swelling rock).

[0180] Each specific component model (such as advanced small duct - high-pressure water-rich fractured zone type) is defined by a set of key parameters, including:

[0181] Geometric parameters: diameter, length, circumferential spacing, longitudinal spacing, and layout angle.

[0182] Material parameters: elastic modulus, yield strength, grouting body permeability, and drainage capacity.

[0183] Behavioral parameters: bond-slip model parameters with surrounding rock, calculation rules for grouting diffusion radius, and drainage efficiency coefficient.

[0184] These key parameters are dynamically adjusted according to the actual engineering requirements (such as adjusting the duct diameter and spacing according to the predicted water inflow).

[0185] For each local risk area identified in step S6, its key features are extracted as the basis for component selection, including:

[0186] Dominant geological anomaly type (DGA): Based on S2 (Geological information integration), S4 (Hierarchical prediction), and S5 (Fracture analysis), the dominant anomaly in the risk area is comprehensively determined. Common types include:

[0187] Fault fracture zone (further distinguish width, filling material, water conductivity), water-rich pocket or high-pressure water inrush area (combining the hydrogeological model of S2 and the fracture connectivity rate of S5), weak interlayer or altered zone, large deformation or squeezing stratum, high in-situ stress rock burst area.

[0188] Surrounding rock stability level (RSL): The probability distribution of the surrounding rock classification in this area output from S4, take the level with the highest probability or the most unfavorable level in the probability distribution of the surrounding rock classification (for example, the probability of grade Ⅳ is 60%, and the probability of grade Ⅴ is 30%, then it is regarded as grade Ⅳ). The higher the level (such as grade Ⅴ), the worse the stability.

[0189] Key quantitative indicators: Such as the predicted maximum water inrush volume, the direction and magnitude of the maximum principal stress (from the initial in-situ stress field), and the fracture connectivity rate threshold (from S5).

[0190] The method of the component selection mapping rule includes: establishing a set of logical decision trees, and mapping the combination of (DGA, RSL, key indicators) to the optimal or recommended component model in the BIM component library;

[0191] Rule examples:

[0192] IF (DGA = F AND RSL = Ⅳ AND fault width > 5m) THEN select steel arch - wide fracture zone strengthened type + pipe shed - long and dense spacing type + full-ring deep-hole grouting body;

[0193] IF (DGA = W AND RSL = Ⅴ AND predicted water inrush volume > 150 m³ / h) THEN select advanced curtain grouting - high-pressure type + drainage hole array - strong drainage type + waterproof board - high-pressure resistant type;

[0194] IF (DGA = D AND RSL = Ⅳ AND predicted convergence > 300mm) THEN select collapsible steel arch + long cable anchor system + temporary invert - heavy type.

[0195] The construction of the BIM component library is based on engineering geological experience, numerical simulation verification and historical case library (the database of S1) to ensure the balance between the safety and economy of the selection; among them, the geological anomaly type (DGA) determines the core coping strategy (such as giving priority to water control, giving priority to preventing collapse), the surrounding rock stability level (RSL) determines the support strength level (such as component size, density, material grade), and the key indicators are used for parameter fine-tuning (such as grouting pressure, bolt length).

[0196] Exemplary description of the intelligent construction process of the local digital twin model:

[0197] In the overall digital twin model of the tunnel, locate the target risk area (such as the section from station number K10+520 to K10+540), accurately delimit its spatial range; determine the dominant geological anomaly type (DGA) for this risk area, such as F (fault fracture zone), determine its surrounding rock stability level (RSL), such as level Ⅳ (probability 65%), extract key indicators, such as fault width = 8m, fracture connectivity rate = 0.85, predicted water inflow = 80m³ / h.

[0198] Input (DGA = F, RSL = Ⅳ, fault width = 8m,...) into the mapping rule engine:

[0199] The engine matches the rules: IF (DGA = F AND RSL = Ⅳ AND fault width > 5m) THEN...

[0200] Output the recommended component combination: [steel arch - enhanced type for wide fracture zone, pipe shed - long and dense - spacing type, full - ring deep - hole grouting body].

[0201] According to the specific dimensions (length, cross - section) of the current risk area and the key indicator values (such as 8m wide), automatically adjust the parameters of the selected components, including:

[0202] Steel arch type: enhanced type for wide fracture zone, parameter adjustment: arch spacing 0.5m (the original library default is 0.75m, encrypted due to large width).

[0203] Pipe shed type: long and dense - spacing type, parameter adjustment: length 30m (covering the fault - affected area), circumferential spacing 0.2m (encrypted due to high connectivity rate).

[0204] Grouting body: full - ring deep - hole grouting body, parameter adjustment: diffusion radius 1.5m (calculated according to fracture aperture and connectivity rate).

[0205] Automatically assemble the parameter - adjusted component models into the three - dimensional spatial position of the risk area. The steel arches are arranged at intervals along the tunnel axis, the pipe sheds are arranged along the excavation contour line according to the designed angle and spacing, and the grouting body model generates equivalent continuous media or discrete grouting units according to the designed range.

[0206] Couple the assembled BIM component system for support - reinforcement - drainage with the refined surrounding rock grid model (including high - resolution geological attributes provided by S4 and S5) in the risk area. The coupling methods include:

[0207] Define the interaction relationship between the components and the surrounding rock (such as seepage - stress coupling between the grouting body and the fracture surface, bond - slip between the bolt and the rock mass);

[0208] Set construction process parameters (such as excavation sequence, support application timing);

[0209] Finally, a parametric BIM local digital twin sub-model that can be directly used for high-precision mechanical-seepage coupling analysis is formed.

[0210] Exemplary:

[0211] Scenario: A high-risk area is identified at station K10+525:

[0212] DGA = W (water-rich pocket): Advanced geophysical prospecting and the S5 fracture model show the presence of a closed high-pressure water body.

[0213] RSL = Grade V: The S4 grading probability field shows that the probability of Grade V at this location reaches 80%.

[0214] Key indicators: Predicted hydrostatic pressure is 2.5 MPa, and water inflow is 200 m³ / h.

[0215] Component selection mapping:

[0216] Rule matching: IF (DGA = W AND RSL = V AND water pressure > 2 MPa) THEN select [advanced curtain grouting - ultra-high pressure type, waterproof board - extra-high pressure resistance type, drainage hole array - extra-large flow rate type].

[0217] Parametric instantiation:

[0218] Advanced curtain grouting - ultra-high pressure type: Automatically set according to the water pressure of 2.5 MPa and fracture characteristics, grouting pressure = 3.0 MPa, and slurry viscosity = low-viscosity chemical slurry.

[0219] Waterproof board - extra-high pressure resistance type: Automatically set thickness = 3 mm, and tensile strength = 30 MPa.

[0220] Drainage hole array - extra-large flow rate type: Automatically set according to the water inflow of 200 m³ / h, hole diameter = 100 mm, hole depth = 6 m, and circumferential spacing = 1.0 m.

[0221] Construction of the local digital twin sub-model: Automatically assemble the above parametric components within the section of K10+520 - 530:

[0222] Generate an ultra-high pressure curtain grouting body (equivalent to a low-permeability reinforcement circle) within a range of 10 m in front of the tunnel face, assemble an extra-high pressure resistance waterproof board inside the primary support layer, and arrange two rows of extra-large flow rate drainage holes at the bottom of the side wall.

[0223] In step S2, it is necessary to integrate multi-source heterogeneous information such as geological drilling data, geophysical exploration data, and construction monitoring data into a unified three-dimensional geological model; achieve sub-millimeter-level spatial positioning and millisecond-level time synchronization through multi-sensor fusion positioning technology to ensure that all data has an accurate spatio-temporal reference.

[0224] Geological borehole data: Spatial positioning relies on a survey control network (accuracy at the centimeter level), and the timestamp is the construction record time (accuracy at the minute level).

[0225] Geophysical exploration data (such as TSP, geological radar): The spatial positions of measurement points are collected in real time by mobile sensors (dynamic positioning is required), the data acquisition frequency is high (such as sampling at the millisecond level), and it is required that the positions of measurement points are strictly synchronized with the detection waveforms.

[0226] Construction monitoring data (such as convergence meters, stress gauges, water inrush sensors): The sensors are fixed on the tunnel wall, and absolute position accuracy (sub - centimeter level) and strict alignment of time scales (millisecond level) for multi - sensor data are required to capture relevant events (such as surrounding rock deformation and stress redistribution after excavation).

[0227] Building a high - precision tunnel digital twin requires all data to achieve sub - millimeter - level spatial accuracy (for fine geological interface characterization and deformation analysis) and millisecond - level time synchronization accuracy (for dynamic process correlation analysis, such as blasting vibration and subsequent deformation) in a unified coordinate system.

[0228] The architecture of multi - sensor fusion positioning technology includes:

[0229] High - precision total station (or laser tracker): Provides static or quasi - static absolute spatial reference (accuracy better than 1mm), used to calibrate the positions of fixed reference stations and key control points (such as tunnel axis points).

[0230] GNSS receiver (Beidou / GPS, etc.): Provides absolute positions (accuracy at the centimeter level) in the global coordinate system and a high - precision UTC time source (nanosecond level) at the tunnel entrance or in open areas on the ground surface, serving as the source of the spatio - temporal reference.

[0231] Inertial measurement unit (IMU): Mounted on mobile detection devices (such as TSP vibration source vehicles, geological radar antenna vehicles), measures acceleration and angular velocity in real time, provides high - frequency (kHz - level) attitude and short - term displacement information, and makes up for the positioning continuity when GNSS fails in the tunnel.

[0232] Ultra - wideband (UWB) indoor positioning system: Fixed base stations are deployed in the tunnel (spacing 50 - 100m), providing relative positioning (accuracy at the centimeter level) for mobile devices (detection devices, handheld terminals) and fixed sensors, and achieving synchronization through timestamp transmission.

[0233] High - precision clock source: Adopts rubidium atomic clocks or industrial switches with PTP (Precision Time Protocol) to provide a unified and stable time reference for all sensors.

[0234] Establish a spatio - temporal reference station at the tunnel entrance, integrating GNSS receivers (acquiring absolute coordinates and UTC), high - precision clock sources, and total stations (calibrating control points at the entrance).

[0235] Deploy a UWB positioning base station network and an optical fiber time synchronization network (transmitting PTP signals) along the tunnel axis.

[0236] Integrate IMU + UWB tags + GNSS (when available at the tunnel entrance) for mobile detection devices.

[0237] Equip fixed monitoring sensors with UWB tags or access time synchronization through a wired network.

[0238] Exemplary description of the method for achieving sub - millimeter - level spatial positioning:

[0239] Use a total station to accurately measure the positions of the tunnel entrance control points and UWB base stations (sub - millimeter level) and incorporate them into the engineering coordinate system.

[0240] The GNSS reference station provides the conversion parameters between the engineering coordinate system and the global coordinate system.

[0241] Dynamic positioning of mobile devices:

[0242] Outside the tunnel or at the tunnel entrance: Integrate GNSS positioning (absolute position) and IMU data (smooth trajectory) with centimeter - level accuracy.

[0243] Inside the tunnel: Switch to UWB positioning (providing the relative distance to the base stations) and integrate IMU data (resolve attitude and velocity), using tightly - coupled Kalman filtering, including:

[0244] State variables: Device position (X, Y, Z), velocity (Vx, Vy, Vz), attitude angles (Roll, Pitch, Yaw), and IMU zero bias.

[0245] Observation variables: UWB ranging values, IMU acceleration / angular velocity.

[0246] Regularly calibrate the UWB base station positions and online calibrate the IMU zero bias through a total station to suppress error accumulation, and finally achieve a dynamic positioning accuracy of <5 mm (sub - millimeter level) for mobile devices inside the tunnel (such as the position of the center point of the ground - penetrating radar antenna).

[0247] Fixed sensor positioning:

[0248] Accurately measure its spatial coordinates (sub - millimeter level) once by a total station during installation.

[0249] Bind its location information to the UWB tag or network address.

[0250] Exemplary (analysis of water inrush events):

[0251] Time T0: The borehole water level gauge located at pile number K10+535.200 (UWB positioning accuracy ±3mm, time synchronization accuracy ±0.5ms) detected a sudden drop in water level (-2.5m).

[0252] Time T0+120ms: The water gushing sensor at K10+520.000, 15m away (positioning accuracy ±2mm, synchronization accuracy ±0.5ms) detected a sudden increase in flow rate of 50m³ / h.

[0253] Spatio-temporal analysis: Combining the accurate spatial distance of 15.200m and time difference of 120ms between the two points, the water flow velocity in the water gushing channel was calculated to be approximately 126.7m / s. Combining this information with the S5 fracture network model can accurately locate the possible dominant water-conducting fractures and provide a target area for water plugging and grouting.

[0254] Obviously, those skilled in the art can make various modifications and variations to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, the present invention is also intended to include these modifications and variations.

[0255] The above are only the preferred specific embodiments of the embodiments of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art within the technical scope disclosed in the present application, according to the technical solution and its concept of the present application, makes equivalent replacements or changes, and all should be covered by the protection scope of the present application.

Claims

1. A spatial registration method for multi-source heterogeneous data in tunnels, characterized in that, Including: Construct a lightweight digital twin model based on the design parameter library, construction log library, and real-time monitoring library; Map the multi-source sensor data in the real-time monitoring library to the spatial nodes of the digital twin model through the spatio-temporal synchronization algorithm; Perform outlier rejection on geological exploration data, drilling parameters, and geological sketches based on the probability distribution model, and use the discrete smooth interpolation algorithm to generate the three-dimensional distribution field of surrounding rock mechanical parameters; Taking the three-dimensional distribution field of surrounding rock mechanical parameters as the input, dynamically predict the probability distribution of the surrounding rock classification ahead through a machine learning model optimized by transfer learning; Based on the three-dimensional point cloud data of the tunnel face in the digital twin model, extract the geometric parameters of the structural plane and construct a fracture network model with normal vector constraints; When the probability of surrounding rock grades I to III in the probability distribution of surrounding rock classification exceeds the set threshold, and the fracture connectivity rate output by the fracture network model exceeds the critical value, associate the support parameters in the design parameter library, the geological event records in the construction log library, and the convergence deformation data to construct a local digital twin sub-model; Perform multi-scale fusion of the local digital twin sub-model and the three-dimensional distribution field of surrounding rock mechanical parameters based on topological similarity to eliminate geometric and attribute discontinuity regions.

2. The spatial registration method for multi-source heterogeneous data in a tunnel according to claim 1, characterized in that For outlier rejection, the Gaussian mixture model is used to verify the consistency of data distribution, and abnormal data deviating from the preset probability interval is rejected.

3. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 1, characterized in that, The machine learning model is a bidirectional LSTM network embedded with a multi-head attention mechanism, and the input parameters include the dynamic coupling characteristics of the rock mass integrity coefficient and the uniaxial compressive strength of the rock.

4. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 1, characterized in that The normal vector constraint specifically includes: identifying the dominant fracture group based on normal vector space clustering, and the spatial autocorrelation characteristics of the fracture density distribution are fused in the clustering process.

5. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 1, characterized in that For multi-scale fusion, the gradient constraint deformation algorithm is used, and the boundary deformation of the local digital twin sub-model is adaptively adjusted according to the spatial gradient of the surrounding rock mechanical parameters in the neighborhood.

6. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 4, characterized in that After identifying the dominant fracture group, use a conditional generative adversarial network to construct a spatial extension probability field of the structural plane and predict the fracture development trend in the unexcavated area ahead.

7. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 5, characterized in that When constructing the local digital twin sub-model, call the parametric BIM component library, and the component selection logic is based on the mapping relationship between the geological anomaly type and the surrounding rock stability level.

8. A spatial registration method for multi-source heterogeneous data in a tunnel according to claim 1, characterized in that The spatio-temporal synchronization algorithm uses multi-sensor fusion positioning technology to achieve sub-millimeter-level spatial positioning and millisecond-level time synchronization.

Images