A tunnel operation area dynamic parameterization modeling method

Abstract

The present application belongs to the technical field of parameterized modeling, and discloses a tunnel operation area dynamic parameterized modeling method; comprising: collecting measured data of a tunnel operation area; based on the collected measured data, establishing an initial model of the tunnel operation area; dividing the initial model of the tunnel operation area into a plurality of model regions, taking the model regions as nodes to construct an operation area undirected graph; assigning colors to the nodes in the operation area undirected graph; for each model region, determining a corresponding parameter group according to the assigned color, and updating to obtain an overall parameterized model based on the initial model of the tunnel operation area and the corresponding parameter groups of the model regions; namely, the dynamic parameterized modeling result of the entire tunnel operation area, which provides a scientific decision basis for design optimization, construction scheme formulation, risk assessment and operation management of tunnel projects.

Description

Technical Field

[0001] This invention relates to the field of parametric modeling technology, and more specifically, to a dynamic parametric modeling method for tunnel operation areas. Background Technology

[0002] Tunnel engineering, as a crucial component of infrastructure construction, plays a vital role in national economic development and social progress. However, due to the complexity and variability of the tunnel operation area environment, existing tunnel operation area modeling techniques suffer from numerous shortcomings, severely restricting the design, construction, and operation management of tunnel projects. Firstly, the multi-source parameters of the tunnel operation area exhibit significant spatial differences, making it difficult to accurately reflect the actual situation using uniform parameters. Taking tunnel engineering in Tibet as an example, the region's complex geological structure, diverse soil and rock properties, and significant impact from the plateau environment necessitate a uniform parameter combination for modeling, inevitably leading to substantial deviations between the model and reality, posing significant risks to engineering design and construction. Secondly, existing tunnel operation area modeling methods are mostly based on simplification assumptions and empirical formulas, lacking full utilization of measured data, resulting in low modeling accuracy and reliability. Designers are forced to rely on limited empirical data, leading to significant discrepancies between the model and reality. Furthermore, traditional modeling methods often treat the tunnel operation area as a whole, failing to employ different parameter combinations for different regions, thus hindering the representation of local characteristics. Finally, improper parameter transition processing can lead to discontinuities or abrupt parameter changes in certain areas, affecting the model's accuracy and reliability.

[0003] In view of this, the present invention proposes a dynamic parametric modeling method for tunnel operation areas to solve the above problems. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art and to achieve the above objectives, the present invention provides the following technical solution: a dynamic parametric modeling method for tunnel operation areas, comprising: step 1, collecting measured data of the tunnel operation area;

[0005] Step 2: Based on the collected measured data, establish an initial model of the tunnel operation area;

[0006] Step 3: Divide the initial model of the tunnel operation area into several model regions, and use the model regions as nodes to construct an undirected graph of the operation area;

[0007] Step 4: Assign colors to the nodes in the undirected graph of the work area;

[0008] Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter group corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

[0009] Furthermore, the measured data of the tunnel operation area includes geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data, and construction parameter data.

[0010] Furthermore, the method for establishing the initial model of the tunnel operation area includes:

[0011] The collected measured data of various types are classified and organized, archived according to data type, and corresponding data packets are obtained. The corresponding data packets are then subjected to smooth interpolation to obtain corresponding smooth data packets.

[0012] The three-dimensional coordinate data of various locations within the tunnel operation area are acquired. The acquired point cloud data is registered and merged into a unified point cloud dataset. The point cloud dataset is then filtered to obtain a standard point cloud dataset. A region-growing-based segmentation algorithm is used to segment the standard point cloud dataset into different point cloud components.

[0013] Fit corresponding surface models to the segmented point cloud data; mesh the fitted surface models to generate a mesh model composed of triangular or quadrilateral patches; optimize the generated mesh model to obtain the preliminary mesh model corresponding to each point cloud data; register and connect the preliminary mesh models corresponding to each point cloud data to obtain the preliminary tunnel operation area model.

[0014] The data points in the smoothing data package are registered to the space of the preliminary tunnel operation area model according to their corresponding coordinate values. In other words, the preliminary tunnel operation area model is given initial parameters to obtain the initial model.

[0015] Furthermore, the method of performing smooth interpolation on the corresponding data packets includes:

[0016] The coordinates of the start and end points of the tunnel line are used as the initial reference points for the spatial range of the data points in the data packet. Based on the initial reference points, the spatial range of the data packet is determined, the number of levels of the spatial range is defined, the entire spatial range is taken as the first level, and the first level is divided into a regular three-dimensional grid.

[0017] Define the order M of the mapping curve as M = ceil(log2(max(n_x,n_y,n_z,α×(-∑(p(i)×log2(p(i)))))))); where n_x, n_y, and n_z are the number of 3D grids in the x, y, and z directions of the 3D space, respectively; p(i) is the probability of the i-th attribute value appearing in the data packet; α is the weight coefficient, and ceil() represents rounding up; based on the defined order M, calculate the encoding range [0, 2] for curve encoding of the data packet. (3×M) -1], meaning the minimum encoded value in the encoding range is 0, and the maximum encoded value is 2. (3×M) -1;

[0018] A virtual mapping curve is constructed based on the order M. The data points in the data packet are mapped onto the mapping curve using the Butz algorithm, and each data point is represented as a curve point. The value of the curve point corresponding to each data point is the encoded value on the mapping curve.

[0019] The encoded values ​​corresponding to data points are divided into several intervals according to the encoding range, and each interval corresponds to a data block. The size of the interval is... Where S is the number of data blocks; the lower and upper bounds of each interval are determined based on the calculated interval size; all data points in the data packet are traversed and assigned to the corresponding intervals according to their encoded values; each interval corresponds to a first-level data block;

[0020] If the number of data points in a data block is greater than the preset data volume threshold, the corresponding data block is treated as a new level and the data is re-divided in the same way as the first level. This process continues until the number of data points in a data block is less than or equal to the preset data volume threshold or the preset number of division levels is reached. Then, the division stops, and the data blocks are merged according to the strategy to obtain the final data block.

[0021] The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations are performed in parallel for the distributed data blocks to obtain the interpolation results of each computing node. The interpolation results from each computing node are then merged to obtain the corresponding smooth data packets.

[0022] Furthermore, the method for merging data blocks according to the strategy includes:

[0023] For each data block, calculate the centroid coordinates of all data points within it; the centroid coordinates represent the spatial center of the data block.

[0024] Construct a spatial index from the centroid coordinates of all data blocks; set an adaptive spatial distance threshold KL = α1 × mean + α2 × std, where α1 and α2 are two distance adjustment coefficients; mean is the average Euclidean distance between the centroid coordinates of all data points within all data blocks; std is the standard deviation of the Euclidean distance between the centroid coordinates of all data points within all data blocks.

[0025] A preset search radius is used to find data blocks within the search radius for each data block, which are then denoted as neighboring data blocks. For each pair of neighboring data blocks, the Euclidean distance between their centroid coordinates is calculated. For two neighboring data blocks whose Euclidean distance is less than the spatial distance threshold, they are merged into a new data block. The merged new data block contains all the data points of the original two data blocks, and the centroid coordinates of the new data block are recalculated. This merging process is repeated until there are no more data block pairs to merge. After each merge, the spatial index and data blocks are updated.

[0026] Furthermore, the fitting method for the surface model includes:

[0027] The point cloud data is divided into a three-dimensional voxel space, and the voxel space is processed by median filtering to obtain the filtered voxel space. Surface voxels are extracted from the filtered voxel space. Surface fitting is performed on the extracted surface voxels to generate a continuous surface model.

[0028] Furthermore, the method for extracting surface voxels from the filtered voxel space includes:

[0029] Define a sliding window to scan the filtered voxel space, calculating the curvature and density features of voxels within each window. Curvature features include principal curvature and Gaussian curvature; density features include voxel occupancy and the number of voxel neighbors. Normalize the curvature and density features and perform weighted calculations to obtain a comprehensive feature value. Based on the distribution of the comprehensive feature value, divide the filtered voxel space into several voxel regions. For each voxel region, define a seed point count Nse. Randomly select Nse seed points for each voxel region.

[0030] Nse = a1 × (Nv) β1 ×exp(γ×σc)×(1+γ1×σd)×(1+δ×σn); where,

[0031] a1, γ, γ1, δ, and β1 are all adjustment coefficients; Nv is the number of voxels in the corresponding voxel region; σc is the standard deviation of the comprehensive eigenvalues ​​in the corresponding voxel region; σd is the standard deviation of the density values ​​of voxels in the corresponding voxel region; and σn is the standard deviation of the normal vectors of voxels in the corresponding voxel region.

[0032] Nse randomly selected seed points are used as the initial region. A curvature threshold is preset. Starting from the randomly selected seed point, the curvature values ​​of its 26 neighboring voxels are calculated. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current initial region and marked as visited. The visited voxel is used as a new seed point, and the curvature values ​​of its 26 neighboring voxels are calculated repeatedly. Based on the curvature values, the voxels to be added to the initial region next are determined. This process continues until no voxels are added to the initial region, resulting in the final initial region. The final initial region is marked as a surface region. All voxels marked as surface regions are collected to form the extracted surface voxels.

[0033] Furthermore, the method of dividing the initial model of the tunnel operation area into several model regions includes:

[0034] The initial model of the tunnel working area was analyzed to identify different structural functional components, including the tunnel body, initial support, permanent lining, pipelines, and ventilation facilities. Each of the initial support, permanent lining, pipelines, and ventilation facilities was considered as a model area.

[0035] The part corresponding to the tunnel body is extracted from the initial model and denoted as the tunnel mesh model. The tunnel mesh model is used as the wave source, and a set of points is uniformly sampled on the surface of the tunnel mesh model as the wave source point set. An initial wavefront is placed at each wave source point, and the simulation time step number is defined. At each simulation time step, the position and shape of the wavefront are updated according to the propagation equation.

[0036] At each time step, the positions and intersection times of the intersection points between different wavefronts are detected; all intersection points are arranged in chronological order of their intersection times, thus obtaining the point set on the centerline; the point set on the centerline is smoothed; the smoothed point set on the centerline is parameterized, thus obtaining the centerline of the tunnel mesh model.

[0037] Spatial sampling is performed along the extracted centerline at preset sampling intervals to obtain several sampling points. Each sampling point corresponds to a cross-sectional location of the tunnel mesh model. At each sampling point, the cross-section of the tunnel mesh model is extracted and denoted as a section. The shape descriptor of each section is extracted. The shape descriptor includes perimeter, area, geometric moment, and central moment. The perimeter, area, geometric moment, and central moment are weighted and summed to obtain the shape feature value corresponding to the section. The burial depth of each section is defined as the distance from the center of the section to the ground surface.

[0038] Initially, each cross section is treated as a separate model region. Preset shape thresholds and burial depth thresholds are used. If the absolute value of the difference between the shape feature values ​​of any two cross sections is less than the shape threshold, and the absolute value of the difference between their burial depth values ​​is less than the burial depth threshold, then the two corresponding cross sections are classified into the same model region.

[0039] Furthermore, methods for constructing an undirected graph of the work area include:

[0040] A unique identifier is assigned to each model region; the geometric center or centroid of each model region is calculated as the representative point of that region; the average Euclidean distance between all pairs of points on the boundary of any two model regions is calculated and denoted as the adjacency metric between the two model regions;

[0041] A pre-defined adjacency threshold is used. If the adjacency metric between two model regions is less than the threshold, the two model regions are considered adjacent. Each model region is treated as a node, and the node's attributes include a unique identifier and the coordinates of the representative point. An empty undirected graph G = (V, E) is constructed, where V is the set of nodes and E is the set of unweighted edges. If two model regions are adjacent, an unweighted edge is added to the undirected graph G, connecting the corresponding nodes of the two model regions. The attributes of the unweighted edge include the adjacency metric between the two model regions. The contact area and intersection volume between two adjacent model regions are calculated as the weight of the unweighted edge. The undirected graph G at this point is the undirected graph of the work area.

[0042] Furthermore, the method for assigning colors to nodes in the undirected graph of the work area includes:

[0043] Calculate the degree of each node in the undirected graph of the work area, define the coloring rule as coloring the node with the highest degree first, and use a greedy strategy to gradually assign color to each node in the undirected graph of the work area until all nodes are assigned color; that is, obtain the initial solution for assigning color to the nodes in the undirected graph of the work area; calculate the number of colors in the initial solution; count the number of node conflicts in the initial solution, that is, the number of adjacent pairs of nodes assigned the same color.

[0044] Define a set of nested neighborhood structures, where each neighborhood structure N_k describes the transition from the current solution to a new solution; k is the index of the neighborhood structure. Randomly select a neighborhood structure from the nested structures and randomly generate a new solution x' within that selected structure. The new solution x' is obtained by locally perturbing the current solution, including swapping node colors and recoloring some nodes. Starting from the new solution x', search for a locally optimal solution x” within the selected neighborhood structure N_k. Calculate the target value of the new solution x', which is the weighted sum of the number of node conflicts and the number of colors. Based on the definition of the selected neighborhood structure, generate a set N(x') of all neighborhood solutions for the new solution x'. Each neighborhood solution in the set N(x') is obtained by locally perturbing x'.

[0045] For each neighborhood solution x”, calculate its objective value and record the currently found optimal neighborhood solution and its objective value; the optimal neighborhood solution is the one with the minimum weighted sum of node conflict number and coloring number; compare the objective value of the new solution x’ with that of the optimal neighborhood solution; if the objective value of the optimal neighborhood solution is less than that of the new solution x’, then move x’ to the optimal neighborhood solution x_best, i.e., x’ = x_best; otherwise, keep x’ unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimum x”;

[0046] If the objective value of the local optimum x” is smaller than the objective value of the initial solution, then update the initial solution to the local optimum x”; otherwise, keep the initial solution unchanged; if no better solution than the local optimum x” is found in the selected neighborhood structure N_k, then replace the neighborhood structure N_k with any of the nested neighborhood structures; and regenerate a new solution, and iterate until the local optimum is not updated for several consecutive times; stop the iteration, output the current optimal solution and its corresponding node coloring scheme; and assign colors to the nodes in the undirected graph of the work area according to the node coloring scheme.

[0047] The technical effects and advantages of the dynamic parametric modeling method for tunnel operation areas proposed in this invention are as follows:

[0048] This invention significantly improves the accuracy, reliability, and applicability of modeling, providing strong technical support for the design, construction, and operation management of tunnel engineering. It fully utilizes comprehensive measured data from the tunnel operation area, including data on geology, mechanics, environment, and construction, providing reliable foundational data for modeling. By cleverly dividing the tunnel operation area into multiple zones and employing different parameter combinations for different zones, it effectively reflects the spatial differences within the tunnel operation area, avoiding errors and deviations caused by using uniform parameters. Furthermore, smooth transition zones are established between adjacent zones, ensuring the continuity and consistency of the overall model. The established dynamic parameterized model not only accurately describes the current state of the tunnel operation area but, more importantly, can dynamically adjust parameters according to actual conditions, predicting and simulating the response of the tunnel operation area under different working conditions. This provides a scientific basis for decision-making in tunnel engineering design optimization, construction plan formulation, risk assessment, and operation management. Attached Figure Description

[0049] Figure 1 This is a schematic diagram of a dynamic parametric modeling method for tunnel operation areas according to the present invention;

[0050] Figure 2 This is a schematic diagram of a dynamic parametric modeling system for tunnel operation areas according to the present invention. Detailed Implementation

[0051] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0052] Example 1

[0053] Please see Figure 1 As shown in this embodiment, a dynamic parametric modeling method for tunnel operation areas includes:

[0054] Step 1: Collect measured data from the tunnel work area;

[0055] Step 2: Based on the collected measured data, establish an initial model of the tunnel operation area;

[0056] Step 3: Divide the initial model of the tunnel operation area into several model regions, and use the model regions as nodes to construct an undirected graph of the operation area;

[0057] Step 4: Assign colors to the nodes in the undirected graph of the work area;

[0058] Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter group corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

[0059] The measured data for the tunnel operation area include geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data, and construction parameter data;

[0060] Geological exploration data includes geological data such as rock and soil strata, lithology, structural planes, and permeability obtained from borehole core sampling and geophysical testing. These data reflect the geological conditions of the tunnel's location. In-situ testing data includes data on rock mechanics parameters and groundwater pressure obtained from in-situ tests such as rock mechanics tests and hydraulic pressure tests. These data reflect the actual mechanical properties of the surrounding rock. Topographic data includes data on topographic relief and elevation along the tunnel route. These data reflect the tunnel's spatial location and burial depth. Stress measurement data includes in-situ stress test data such as horizontal stress and vertical stress. These data reflect the stress environment of the tunnel. Environmental parameter data includes environmental parameter test data such as temperature and humidity. These data affect the mechanical properties of the surrounding rock. Construction parameter data includes construction parameter data such as excavation methods and support measures. These data are important boundary conditions for establishing the initial model.

[0061] The above measured data comprehensively reflects the geological, mechanical, environmental and construction conditions of the tunnel operation area, and is the basic data for establishing the initial model of the tunnel operation area.

[0062] The methods for establishing an initial model of the tunnel work area include:

[0063] The collected measured data of various types are classified and organized, and archived according to data type (such as geological data, mechanical data, environmental data, etc.) to obtain corresponding data packets. The corresponding data packets are then subjected to smooth interpolation to obtain the corresponding continuous data distribution; that is, the corresponding smooth data packets.

[0064] Specifically, the coordinates of the tunnel's origin and destination points are used as preliminary reference points for the spatial range of data points within the data packet. Based on these preliminary reference points, the spatial range of the data packet is determined. The spatial range is defined as the range formed by the minimum and maximum coordinate values ​​of the three coordinate systems in three-dimensional space. The number of levels for dividing the spatial range is defined (a suitable number of levels is determined through trial and error). The entire spatial range is taken as the first level, and the first level is divided into a regular three-dimensional grid (the number of three-dimensional grids is determined based on the preset number of data blocks). The order of the mapping curve is defined as M = ceil(log2(max(n_x,n_y)). ,n_z,α×(-∑(p(i)×log2(p(i)))))))); where n_x, n_y, and n_z are the number of three-dimensional grids in the x, y, and z directions of the three-dimensional space, respectively; p(i) is the probability of the i-th attribute value appearing in the data packet, calculated by dividing the number of times the attribute value appears by the number of data points in the data packet; the attribute value is any one of the x, y, and z coordinates of the data point; α is the weight coefficient, which is adjusted according to the actual situation, and ceil() means rounding up; according to the defined order M, the encoding range [0, 2] of the data packet for curve encoding is calculated. (3×M) -1], meaning the minimum encoded value in the encoding range is 0, and the maximum encoded value is 2. (3×M) -1.

[0065] A virtual mapping curve is constructed based on the order M. Data points within the data packet are mapped onto the mapping curve using the Butz algorithm, and each data point is represented as a curve point. The value of the curve point corresponding to each data point is the encoded value on the mapping curve.

[0066] The encoded values ​​corresponding to data points are divided into several intervals according to the encoding range, and each interval corresponds to a data block. The size of the interval is... Where S is the number of data blocks; based on the calculated size of the interval, determine the lower and upper bounds of each interval; specifically, the lower bound of the first interval is the minimum encoded value, and the upper bound is the minimum encoded value + the size of the interval - 1. The lower bound of subsequent intervals is the upper bound of the previous interval + 1, and the upper bound is the lower bound + the size of the interval - 1.

[0067] Traverse all data points within the data packet and assign them to corresponding intervals based on their encoded values. Each interval corresponds to a first-level data block, containing all data points falling within that interval. Obtain the number of data points within the data block. If it exceeds a preset data volume threshold, treat the corresponding data block as a new level and re-divide the data in the same way as the first level. Continue until the number of data points within the data block is less than or equal to the preset data volume threshold or the preset number of division levels is reached. Stop the division and merge the data blocks according to the strategy to obtain the final data block.

[0068] Specifically, for each data block, the centroid coordinates of all data points within it are calculated (approximate using the mean of the data point coordinates); the centroid coordinates represent the spatial center of the data block; a spatial index, such as an R-tree or quadtree, is constructed from the centroid coordinates of all data blocks; an adaptive spatial distance threshold KL = α1 × mean + α2 × std is set, where α1 and α2 are two distance adjustment coefficients used to control the size of the threshold; mean is the average Euclidean distance between the centroid coordinates of all data points within all data blocks; and std is the standard deviation of the Euclidean distance between the centroid coordinates of all data points within all data blocks.

[0069] A preset search radius is used to find data blocks within the search radius for each data block, which are then denoted as neighboring data blocks. For each pair of neighboring data blocks, the Euclidean distance between their centroid coordinates is calculated. For two neighboring data blocks whose Euclidean distance is less than the spatial distance threshold, they are merged into a new data block. The merged new data block contains all the data points of the original two data blocks, and the centroid coordinates of the new data block are recalculated. This merging process is repeated until there are no more data block pairs to merge. After each merge, the spatial index and data blocks are updated.

[0070] The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations, such as inverse distance weighted interpolation, kriging interpolation, or radial basis function interpolation, are performed in parallel for the distributed data blocks. The interpolation results of each computing node are obtained. The interpolation results of each computing node are merged to complete the smooth interpolation process, and the corresponding smooth data packets are obtained.

[0071] The three-dimensional coordinate data of each location point in the tunnel operation area are acquired, the acquired point cloud data are registered and merged into a unified point cloud dataset, the point cloud dataset is filtered to obtain a standard point cloud dataset, and a region growing-based segmentation algorithm is used to segment the standard point cloud dataset into different point cloud parts.

[0072] The segmented point cloud data is fitted with corresponding surface models; the fitted surface models are then meshed to generate a mesh model composed of triangular or quadrilateral patches; the generated mesh models are then smoothed and simplified to improve model quality and computational efficiency; a preliminary mesh model corresponding to each point cloud data segment is obtained; the preliminary mesh models corresponding to each point cloud data segment are registered and connected to obtain a preliminary tunnel operation area model; the data points in the smoothed data package are registered to the space where the preliminary tunnel operation area model is located according to their corresponding coordinate values, that is, the preliminary tunnel operation area model is assigned initial parameters; the initial model is obtained.

[0073] The fitting methods for surface models include:

[0074] The point cloud data is divided into a three-dimensional voxel space, and each voxel records whether the spatial location is occupied by the point cloud data. The voxel space is then subjected to median filtering to remove isolated voxels and noisy voxels, resulting in the filtered voxel space.

[0075] Surface voxels are extracted from the filtered voxel space. Specifically, a sliding window is defined, and the filtered voxel space is scanned using this sliding window. The curvature and density features of voxels within each sliding window are calculated. The curvature features include principal curvature and Gaussian curvature. The density features include voxel occupancy (how many voxels are occupied or filled by point cloud data in a given 3D voxel grid) and voxel neighborhood (how many occupied neighbor voxels surround each occupied voxel in a given voxel grid). The curvature and density features are then normalized. The filtered voxel space is divided into several voxel regions based on the distribution of the comprehensive feature values. Specifically, the K-Means algorithm is used to divide the voxels into K clusters based on their comprehensive feature values. The mean difference of the comprehensive feature values ​​between any two clusters is calculated, and a cluster distance threshold is set. If the mean difference between two clusters is less than the cluster distance threshold, they are merged into a new cluster. This merging process is repeated until there are no more cluster pairs to merge, and the final cluster is the voxel region. Each voxel region represents a different combination of curvature and density features.

[0076] For each voxel region, define a seed point count Nse; for each voxel region, randomly select Nse seed points.

[0077] Nse = a1 × (Nv) β1 ×exp(γ×σc)×(1+γ1×σd)×(1+δ×σn); where,

[0078] a1, γ, γ1, δ, and β1 are all adjustment coefficients; Nv is the number of voxels in the corresponding voxel region; σc is the standard deviation of the comprehensive eigenvalues ​​in the corresponding voxel region; σd is the standard deviation of the density values ​​of voxels in the corresponding voxel region, reflecting the degree of density change; σn is the standard deviation of the normal vectors of voxels in the corresponding voxel region, reflecting the degree of normal vector change.

[0079] Nse randomly selected seed points are used as the initial region. A curvature threshold is preset. Starting from the randomly selected seed point, the curvature values ​​of its 26 neighboring voxels are calculated. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current initial region and marked as visited. The visited voxel is used as a new seed point, and the curvature values ​​of its 26 neighboring voxels are calculated repeatedly. Based on the curvature values, the next voxels to be added to the initial region are determined. This process continues until no voxels are added to the initial region, resulting in the final initial region. The final initial region is marked as a surface region. All voxels marked as surface regions are collected to form the extracted surface voxels. Surface fitting is performed on the extracted surface voxels to generate a continuous surface model. Surface fitting methods include moving least squares and Gaussian kernel regression.

[0080] The methods for dividing the initial model of the tunnel operation area into several model regions include:

[0081] The initial model of the tunnel working area was analyzed to identify different structural functional components, including the tunnel body, initial support, permanent lining, pipelines, and ventilation facilities. Each of the initial support, permanent lining, pipelines, and ventilation facilities was considered as a model area.

[0082] The part corresponding to the tunnel body is extracted from the initial model and denoted as the tunnel mesh model. The tunnel mesh model is used as the wave source, and a set of points is uniformly sampled on the surface of the tunnel mesh model as the wave source point set. An initial wavefront is placed at each wave source point. The initial wavefront can be spherical, planar, or other simple shapes. The simulation time step is defined. At each simulation time step, the position and shape of the wavefront are updated according to the propagation equation. The propagation equation is a cascade equation or a fast-running sub-equation.

[0083] At each time step, the positions and intersection times of intersection points between different wavefronts are detected; all intersection points are arranged in chronological order of intersection time, thus obtaining the point set on the centerline; the point set on the centerline is smoothed (e.g., by moving average, B-Slide curve fitting, etc.); the smoothed point set on the centerline is parametrically represented (e.g., by B-Slide curve, piecewise linear curve); thus obtaining the centerline of the tunnel mesh model.

[0084] Spatial sampling is performed along the extracted centerline at preset sampling intervals to obtain several sampling points. Each sampling point corresponds to a cross-sectional position of the tunnel mesh model. At each sampling point, the cross-section of the tunnel mesh model is extracted (using methods such as plane cutting and ray projection for cross-sectional extraction), and denoted as the cross-section. The shape descriptor of each cross-section is extracted. The shape descriptor includes perimeter, area, geometric moments, and central moments. The perimeter, area, geometric moments, and central moments are weighted and summed to obtain the shape feature value corresponding to the cross-section. The burial depth of each cross-section is defined as the distance from the center of the cross-section to the ground surface. The geometric moment refers to a set of values ​​obtained by integrating the coordinate values ​​of each point in a region after multiplying them by the corresponding weights. The central moment is a set of moment values ​​calculated by first translating the origin of the shape's coordinates to the centroid of the shape based on the calculation of the geometric moments.

[0085] Initially, each cross section is treated as a separate model region. Preset shape thresholds and burial depth thresholds are used. If the absolute value of the difference between the shape feature values ​​of any two cross sections is less than the shape threshold, and the absolute value of the difference between their burial depth values ​​is less than the burial depth threshold, then the two corresponding cross sections are classified into the same model region.

[0086] Assign a unique identifier (ID) to each model region; calculate the geometric center or centroid of each model region as the representative point of the region; calculate the average Euclidean distance between all pairs of points on the boundary of any two model regions; record this as the adjacency metric between the two model regions; preset an adjacency threshold, and if the adjacency metric between two model regions is less than the adjacency threshold, then the two model regions are adjacent.

[0087] Each model region is treated as a node, with attributes including a unique identifier and coordinates representing the point. An empty undirected graph G = (V, E) is constructed, where V is the set of nodes and E is the set of unweighted edges. If two model regions are adjacent, an unweighted edge is added to the undirected graph G, connecting the corresponding nodes of the two model regions. The attributes of the unweighted edge include the adjacency metric between the two model regions. The contact area and intersection volume between two adjacent model regions are calculated as the weight of the unweighted edge. The undirected graph G at this point is the job area undirected graph.

[0088] The ways to assign colors to nodes in an undirected graph of a work area include:

[0089] Calculate the degree of each node in the undirected graph of the work area, define the coloring rule as prioritizing the coloring of the node with the highest degree, and use a greedy strategy to gradually assign color to each node in the undirected graph of the work area until all nodes are assigned color; thus, we obtain the initial solution (node ​​coloring scheme) for assigning color to the nodes in the undirected graph of the work area; calculate the number of colors used in the initial solution, i.e., the number of color types used; count the number of node conflicts in the initial solution, i.e., the number of adjacent pairs of nodes assigned the same color.

[0090] Define a set of nested neighborhood structures, where each neighborhood structure N_k describes the transition from the current solution to the new solution; k is the index of the neighborhood structure; for example, N_1 swaps the colors of two adjacent nodes; N_2 recolors a node and all its adjacent nodes; N_3 recolors all nodes in a maximal clique; and N_k_max generates a completely new coloring scheme.

[0091] In a set of nested neighborhood structures, a neighborhood structure is randomly selected, and a new solution x' is randomly generated in the selected neighborhood structure. The new solution x' is obtained by locally perturbing the current solution. The local perturbation methods include swapping node colors and recoloring some nodes.

[0092] Starting from the new solution x', a local optimum x” is found in the selected neighborhood structure N_k; specifically, the target value of the new solution x' is calculated, which is the weighted sum of the number of node conflicts and the number of colors. According to the definition of the selected neighborhood structure, the set of all neighborhood solutions N(x') of the new solution x' is generated; each neighborhood solution in the set of neighborhood solutions N(x') is obtained by locally perturbing x'.

[0093] For each neighborhood solution x", calculate its objective value (the weighted sum of the number of node conflicts and the number of colors), and record the currently found optimal neighborhood solution and its objective value; the optimal neighborhood solution is the one with the minimum weighted sum of the number of node conflicts and the number of colors.

[0094] Compare the objective value of the new solution x' with that of the optimal neighborhood solution; if the objective value of the optimal neighborhood solution is less than that of the new solution x', then move x' to the optimal neighborhood solution x_best, i.e., x' = x_best; otherwise, keep x' unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimum x.

[0095] If the objective value of the local optimum x” is smaller than the objective value of the initial solution, then update the initial solution to the local optimum x”; otherwise, keep the initial solution unchanged; if no better solution (smaller objective value) than the local optimum x” is found in the selected neighborhood structure N_k, then replace the neighborhood structure N_k with any of the nested neighborhood structures; and regenerate a new solution and iterate until the local optimum is not updated for several consecutive times; stop the iteration, and output the current optimal solution (minimum objective value) and its corresponding node coloring scheme.

[0096] According to the node coloring scheme, the nodes in the undirected graph of the work area are assigned colors. In the above process, the neighborhood structure determines the search range and accuracy, and it is usually necessary to gradually expand from a small range to a large range. By systematically exploring different neighborhoods, it is possible to effectively escape local optima and obtain a better node coloring and parameter allocation scheme.

[0097] Based on the color assigned to each model region, several different parameter sets are prepared in advance, each containing several specific parameters (adjusted based on measured data from the tunnel operation area) to comprehensively describe the geological, mechanical, environmental, and construction conditions of the corresponding region. All model regions are traversed, and for the j-th model region, its assigned color c_j is obtained. Among the several different parameter sets, the parameter set corresponding to color c_j is found, and this parameter set is determined as the parameter set corresponding to the j-th model region.

[0098] At the boundary between two adjacent model regions, an overlapping region is set. Within the overlapping region, interpolation is performed on the parameter sets of the two adjacent model regions to achieve a smooth transition. The interpolation method can be linear interpolation, nonlinear interpolation (such as Gaussian process regression), etc.

[0099] Through the above steps, each model area in the entire tunnel operation area obtained a corresponding parametric model, and the parameters were smoothly transitioned at the boundaries of the areas, thus completing the dynamic parametric modeling of the entire operation area.

[0100] The initial model is split according to the division of model regions to obtain initial sub-models corresponding to each model region. For each initial sub-model of the model region, the relevant parameters in the initial sub-model are updated according to the parameter group corresponding to the model region to obtain the parameterized sub-model. For example, the geological model parameters, such as soil and rock strata, lithology, and structural planes, are updated; mechanical parameters, such as rock strength and deformation modulus, are updated; environmental parameters, such as temperature and humidity, are updated; and construction parameters, such as excavation method and support measures, are updated.

[0101] The parametric sub-model is concatenated with the initial model to obtain a global parametric model that includes all model regions. For parametric models in adjacent regions, the overlapping region method is used for smooth interpolation to achieve seamless parameter transition. It is important to note that during parameter updates, the correlation and constraints between parameters must be fully considered to ensure the rationality and consistency of the parameters. At the same time, the parameters need to be corrected and optimized according to the actual situation to improve the accuracy and reliability of modeling.

[0102] This embodiment significantly improves the accuracy, reliability, and applicability of modeling, providing strong technical support for the design, construction, and operation management of tunnel engineering. It fully utilizes comprehensive measured data from the tunnel operation area, including data on geology, mechanics, environment, and construction, providing reliable foundational data for modeling. By cleverly dividing the tunnel operation area into multiple regions and using different parameter combinations for different regions, it effectively reflects the spatial differences within the tunnel operation area, avoiding errors and deviations caused by using uniform parameters. Simultaneously, smooth transition zones are set between adjacent regions, ensuring the continuity and consistency of the overall model. The established dynamic parameterized model not only accurately describes the current state of the tunnel operation area but, more importantly, can dynamically adjust parameters according to actual conditions, predicting and simulating the response of the tunnel operation area under different working conditions. This provides a scientific basis for decision-making in tunnel engineering design optimization, construction plan formulation, risk assessment, and operation management.

[0103] Example 2

[0104] Please see Figure 2 As shown, parts not described in detail in this embodiment are described in Embodiment 1. A dynamic parametric modeling system for tunnel operation areas is provided, including:

[0105] The data acquisition module is used to collect measured data from the tunnel working area.

[0106] The preliminary modeling module establishes an initial model of the tunnel operation area based on the collected measured data;

[0107] The region division and graphing module is used to divide the initial model of the tunnel operation area into several model regions, and use the model regions as nodes to construct an undirected graph of the operation area.

[0108] The coloring module is used to assign colors to nodes in the undirected graph of the work area;

[0109] The parameterization module is used to determine the corresponding parameter group for each model area based on its assigned color. Based on the initial model of the tunnel operation area and the parameter groups corresponding to each model area, the overall parameterized model is updated to obtain the dynamic parameterized modeling result of the entire tunnel operation area. The modules are connected to each other via wired and / or wireless means to realize data transmission between modules.

[0110] Example 3

[0111] This embodiment discloses an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the operation mode of the dynamic parameterized modeling method for tunnel operation areas described above.

[0112] Since the electronic device described in this embodiment is the electronic device used to implement the dynamic parametric modeling method for tunnel operation areas described in this application, those skilled in the art can understand the specific implementation method and various variations of the electronic device in this embodiment based on the dynamic parametric modeling method for tunnel operation areas described in this application. Therefore, how the electronic device implements the method in this application will not be described in detail here. Any electronic device used by those skilled in the art to implement the dynamic parametric modeling method for tunnel operation areas described in this application falls within the scope of protection of this application.

[0113] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters and thresholds in the formulas are set by those skilled in the art according to the actual situation.

[0114] The above description is merely a preferred embodiment of the present invention, and the scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for users of ordinary technical skills, any improvements and modifications made without departing from the principles of the present invention should also be considered within the scope of protection of the present invention.

Claims

1. A method of dynamic parametric modeling of a tunneling work zone, characterized in that, include: Step 1: Collect measured data from the tunnel work area; Step 2: Based on the collected measured data, establish an initial model of the tunnel operation area; Step 3: Divide the initial model of the tunnel operation area into several model regions, and use the model regions as nodes to construct an undirected graph of the operation area; The method of dividing the initial model of the tunnel operation area into several model regions includes: The initial model of the tunnel working area was analyzed to identify different structural functional components, including the tunnel body, initial support, permanent lining, pipelines, and ventilation facilities. Each of the initial support, permanent lining, pipelines, and ventilation facilities was considered as a model area. The part corresponding to the tunnel body is extracted from the initial model and denoted as the tunnel mesh model. The tunnel mesh model is used as the wave source, and a set of points is uniformly sampled on the surface of the tunnel mesh model as the wave source point set. An initial wavefront is placed at each wave source point, and the simulation time step number is defined. At each simulation time step, the position and shape of the wavefront are updated according to the propagation equation. At each time step, the positions and intersection times of the intersection points between different wavefronts are detected; all intersection points are arranged in chronological order of their intersection times, thus obtaining the point set on the centerline; the point set on the centerline is smoothed; the smoothed point set on the centerline is parameterized, thus obtaining the centerline of the tunnel mesh model. Spatial sampling is performed along the extracted centerline at preset sampling intervals to obtain several sampling points. Each sampling point corresponds to a cross-sectional location of the tunnel mesh model. At each sampling point, the cross-section of the tunnel mesh model is extracted and denoted as a section. The shape descriptor of each section is extracted. The shape descriptor includes perimeter, area, geometric moment, and central moment. The perimeter, area, geometric moment, and central moment are weighted and summed to obtain the shape feature value corresponding to the section. The burial depth of each section is defined as the distance from the center of the section to the ground surface. Initially, each cross section is treated as a separate model region. Preset shape threshold and burial depth threshold. If the absolute value of the difference between the shape feature values ​​of any two cross sections is less than the shape threshold, and the absolute value of the difference between their burial depth values ​​is less than the burial depth threshold, then the corresponding two cross sections are classified into the same model region. Methods for constructing an undirected graph of the work area include: A unique identifier is assigned to each model region; the geometric center or centroid of each model region is calculated as the representative point of that region; the average Euclidean distance between all pairs of points on the boundary of any two model regions is calculated and denoted as the adjacency metric between the two model regions; A pre-defined adjacency threshold is used. If the adjacency metric between two model regions is less than the threshold, the two model regions are considered adjacent. Each model region is treated as a node, and the node's attributes include a unique identifier and the coordinates of the representative point. An empty undirected graph G=(V,E) is constructed, where V is the set of nodes and E is the set of unweighted edges. If two model regions are adjacent, an unweighted edge is added to the undirected graph G, connecting the corresponding nodes of the two model regions. The attributes of the unweighted edge include the adjacency metric between the two model regions. The contact area and intersection volume between two adjacent model regions are calculated as the weight of the unweighted edge. The undirected graph G at this point is the undirected graph of the work area. Step 4: Assign colors to the nodes in the undirected graph of the work area; Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter group corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

2. The method of claim 1, wherein, The measured data for the tunnel operation area include geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data, and construction parameter data.

3. The method of claim 2, wherein, The methods for establishing the initial model of the tunnel operation area include: The collected measured data of various types are classified and organized, archived according to data type, and corresponding data packets are obtained. The corresponding data packets are then subjected to smooth interpolation to obtain corresponding smooth data packets. The three-dimensional coordinate data of various locations within the tunnel operation area are acquired. The acquired point cloud data is registered and merged into a unified point cloud dataset. The point cloud dataset is then filtered to obtain a standard point cloud dataset. A region-growing-based segmentation algorithm is used to segment the standard point cloud dataset into different point cloud components. Fit corresponding surface models to the segmented point cloud data; mesh the fitted surface models to generate a mesh model composed of triangular or quadrilateral patches; optimize the generated mesh model to obtain the preliminary mesh model corresponding to each point cloud data; register and connect the preliminary mesh models corresponding to each point cloud data to obtain the preliminary tunnel operation area model. The data points in the smoothing data package are registered to the space of the preliminary tunnel operation area model according to their corresponding coordinate values. In other words, the preliminary tunnel operation area model is given initial parameters to obtain the initial model.

4. The method of claim 3, wherein, The method of performing smooth interpolation on the corresponding data packets includes: The coordinates of the start and end points of the tunnel line are used as the initial reference points for the spatial range of the data points in the data packet. Based on the initial reference points, the spatial range of the data packet is determined, the number of levels of the spatial range is defined, the entire spatial range is taken as the first level, and the first level is divided into a regular three-dimensional grid. Define the order of the mapping curve ;in, , and These represent the number of 3D grids in the x, y, and z directions of the 3D space, respectively. For the first The probability of each attribute value appearing in the data packet; These are the weighting coefficients. Indicates rounding up; based on the defined order. Calculate the encoding range [0, ..., for curve encoding of the data packet] ], meaning the minimum encoded value is 0 and the maximum encoded value is 0. ; Based on order Construct a virtual mapping curve and map the data points in the data packet onto the mapping curve using the Butz algorithm, which is a curve point. The value of the curve point corresponding to each data point is the encoded value on the mapping curve. The encoded values ​​corresponding to data points are divided into several intervals according to the encoding range, and each interval corresponds to a data block. The size of the interval is... ,in The number of data blocks is determined; based on the calculated size of the interval, the lower and upper bounds of each interval are determined; all data points within the data packet are traversed and assigned to the corresponding intervals according to their encoded values; each interval corresponds to a first-level data block; If the number of data points in a data block is greater than the preset data volume threshold, the corresponding data block is treated as a new level and the data is re-divided in the same way as the first level. This process continues until the number of data points in a data block is less than or equal to the preset data volume threshold or the preset number of division levels is reached. Then, the division stops, and the data blocks are merged according to the strategy to obtain the final data block. The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations are performed in parallel for the distributed data blocks to obtain the interpolation results of each computing node. The interpolation results from each computing node are then merged to obtain the corresponding smooth data packets.

5. The method for dynamic parametric modeling of tunnel operation areas according to claim 4, characterized in that, The methods for merging data blocks according to the strategy include: For each data block, calculate the centroid coordinates of all data points within it; the centroid coordinates represent the spatial center of the data block. Construct a spatial index from the centroid coordinates of all data blocks; set an adaptive spatial distance threshold. ,in, and These are two distance adjustment coefficients; The average Euclidean distance between the centroid coordinates of all data points within all data blocks; is the standard deviation of the Euclidean distance between the centroid coordinates of all data points within all data blocks; A preset search radius is used to find data blocks within the search radius for each data block, which are then denoted as neighboring data blocks. For each pair of neighboring data blocks, the Euclidean distance between their centroid coordinates is calculated. For two neighboring data blocks whose Euclidean distance is less than the spatial distance threshold, they are merged into a new data block. The merged new data block contains all the data points of the original two data blocks, and the centroid coordinates of the new data block are recalculated. This merging process is repeated until there are no more data block pairs to merge. After each merge, the spatial index and data blocks are updated.

6. The method for dynamic parametric modeling of tunnel operation areas according to claim 5, characterized in that, The fitting methods for the surface model include: The point cloud data is divided into a three-dimensional voxel space, and the voxel space is processed by median filtering to obtain the filtered voxel space. Surface voxels are extracted from the filtered voxel space. Surface fitting is performed on the extracted surface voxels to generate a continuous surface model.

7. The method for dynamic parametric modeling of tunnel working areas according to claim 6, characterized in that, The method for extracting surface voxels from the filtered voxel space includes: Define a sliding window to scan the filtered voxel space, calculating the curvature and density features of voxels within each window. Curvature features include principal curvature and Gaussian curvature; density features include voxel occupancy and the number of voxel neighbors. Normalize the curvature and density features and perform weighted calculations to obtain a comprehensive feature value. Based on the distribution of the comprehensive feature value, divide the filtered voxel space into several voxel regions. For each voxel region, define a seed point count. For each voxel region, select randomly. Seed points; ;in, , , , and All are adjustment coefficients; This represents the number of voxels within the corresponding voxel region. This represents the standard deviation of the comprehensive eigenvalues ​​within the corresponding voxel region; This represents the standard deviation of the density values ​​of voxels within the corresponding voxel region; This represents the standard deviation of the normal vector of the voxel within the corresponding voxel region; Randomly selected Using seed points as initial regions, and with a preset curvature threshold, the curvature values ​​of the voxels in the 26 neighborhoods of a randomly selected seed point are calculated. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current initial region and marked as visited. The visited voxel is then used as a new seed point, and the curvature values ​​of its 26 neighborhood voxels are calculated repeatedly. Based on the curvature values, the voxels to be added to the initial region are determined. This process continues until no voxels are added to the initial region, resulting in the final initial region. This final initial region is then marked as a surface region. All voxels marked as surface regions are collected together to form the extracted surface voxels.

8. The method for dynamic parametric modeling of tunnel working areas according to claim 7, characterized in that, The methods for assigning colors to nodes in the undirected graph of the work area include: Calculate the degree of each node in the undirected graph of the work area, define the coloring rule as coloring the node with the highest degree first, and use a greedy strategy to gradually assign color to each node in the undirected graph of the work area until all nodes are assigned color; that is, obtain the initial solution for assigning color to the nodes in the undirected graph of the work area; calculate the number of colors in the initial solution; count the number of node conflicts in the initial solution, that is, the number of adjacent pairs of nodes assigned the same color. Define a set of nested neighborhood structures, where each neighborhood structure N_k describes the transition from the current solution to a new solution; k is the index of the neighborhood structure. Randomly select a neighborhood structure from the nested structures and randomly generate a new solution x' within that structure. The new solution x' is obtained by locally perturbing the current solution, including swapping node colors and recoloring some nodes. Starting from the new solution x', search for a locally optimal solution x'' within the selected neighborhood structure N_k. Calculate the target value of the new solution x', which is the weighted sum of the number of node conflicts and the number of colors. Based on the definition of the selected neighborhood structure, generate a set N(x') of all neighborhood solutions for the new solution x'. Each neighborhood solution in the set N(x') is obtained by locally perturbing x'. For each neighborhood solution x'', calculate its objective value and record the currently found optimal neighborhood solution and its objective value; the optimal neighborhood solution is the one with the minimum weighted sum of node conflict number and coloring number; compare the objective value of the new solution x' with that of the optimal neighborhood solution; if the objective value of the optimal neighborhood solution is less than that of the new solution x', then move x' to the optimal neighborhood solution x_best, i.e., x' = x_best; otherwise, keep x' unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimum x''; If the objective value of the local optimum x'' is smaller than the objective value of the initial solution, then update the initial solution to the local optimum x''; otherwise, keep the initial solution unchanged. If no better solution than the local optimum x'' is found in the selected neighborhood structure N_k, then replace the neighborhood structure N_k with any of the nested neighborhood structures; and regenerate a new solution, and iterate until the local optimum is not updated for several consecutive times; stop the iteration, output the current optimal solution and its corresponding node coloring scheme; and assign colors to the nodes in the undirected graph of the work area according to the node coloring scheme.

Images