Three-dimensional dynamic simulation method and system for whole process of tunnel blasting

Abstract

The present application relates to the technical field of three-dimensional simulation, and particularly relates to a tunnel blasting whole-process three-dimensional dynamic simulation method and system. The method comprises the following steps: obtaining a tunnel three-dimensional design digital model; identifying a tunnel high-risk geometric region in the tunnel three-dimensional design digital model; constructing a space blasting constraint field covering a blasting rock mass based on the tunnel high-risk geometric region, and planning an ideal blasting borehole path; deducing a blasting energy sequence according to the ideal blasting borehole path, and performing rock mass damage dynamics blasting simulation with the blasting energy sequence as a boundary condition to evaluate a borehole layout correction parameter for correcting underbreak risk; and applying the borehole layout correction parameter to a preset tunnel blasting scheme to obtain an optimized tunnel blasting scheme. The present application realizes accurate dynamic simulation of a non-circular tunnel blasting process by constructing a space blasting constraint field, so as to solve the problem of inaccurate prediction of overbreak or underbreak at key positions of a complex section.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional simulation technology, and in particular to a three-dimensional dynamic simulation method and system for the entire process of tunnel blasting. Background Technology

[0002] In projects such as high-speed railway tunnels and urban subway stations, tunnel cross-sections are often no longer standard circles or rectangles, but rather employ more complex non-circular profiles, such as composite forms composed of various combined curves, including three-centered circles, five-centered circles, curved sidewalls, straight walls, and arches. Ideal perimeter hole blasting (also known as smooth blasting) requires that the blast holes closely conform to the design profile to form a flat excavation boundary and minimize disturbance to the surrounding rock. However, in non-circular profiles, the radius of curvature of the profile continuously changes from the arch crown and waist to the sidewalls and arch feet, presenting a complex dynamic relationship in three-dimensional space. However, existing simulation methods generally simplify complex curved contours into standard circular or rectangular geometries. While this simplification reduces the difficulty of modeling, it comes at the cost of losing crucial spatial geometric information on the contour, especially the true shape of locations with drastic curvature changes, such as the arch waist and arch foot. Along the simplified regular contour lines, the peripheral holes are evenly arranged at equal intervals or angles, completely ignoring the inherent constraints of contour curvature changes on the density of borehole arrangement, resulting in blasting schemes that are out of touch with actual engineering needs. Summary of the Invention

[0003] Based on this, the present invention provides a three-dimensional dynamic simulation method and system for the entire process of tunnel blasting, in order to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a three-dimensional dynamic simulation method for the entire tunnel blasting process is provided, comprising the following steps: Step S1: Obtain a 3D design digital model of the tunnel with a non-circular contour; identify high-risk geometric areas of the tunnel in the 3D design digital model; Step S2: Construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric area of ​​the tunnel; plan the ideal blasting borehole path by analyzing the cumulative weight of the blasting constraint values ​​in the spatial blasting constraint field; determine the key points of simulation monitoring according to the high-risk geometric area of ​​the tunnel. Step S3: Obtain the three-dimensional geological model of the tunnel; derive the blasting energy sequence based on the ideal blasting borehole path, and use this as the boundary condition to perform rock mass damage dynamics blasting simulation on the three-dimensional geological model of the tunnel. By tracking the energy transfer attenuation rate of key monitoring points in the simulation, evaluate the borehole layout correction parameters used to correct the risk of under-excavation. Step S4: Apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the optimization of the blasting scheme and obtain the optimized tunnel blasting scheme.

[0005] This invention also provides a three-dimensional dynamic simulation system for the entire tunnel blasting process, which executes the three-dimensional dynamic simulation method for the entire tunnel blasting process as described above. The three-dimensional dynamic simulation system for the entire tunnel blasting process includes: The geometric risk identification module is used to acquire a 3D design model of a tunnel with a non-circular contour; and to identify high-risk geometric areas in the 3D design model of the tunnel. The borehole path planning module is used to construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric area of ​​the tunnel; by analyzing the cumulative weight of the blasting constraint values ​​in the spatial blasting constraint field, it plans the ideal blasting borehole path; and determines the key points of simulation monitoring based on the high-risk geometric area of ​​the tunnel. The blasting effect simulation module is used to acquire a three-dimensional geological model of the tunnel; blasting energy sequence is derived based on the ideal blasting borehole path, and rock mass damage dynamics blasting simulation is performed on the three-dimensional geological model of the tunnel using this as boundary condition; by tracking the energy transfer attenuation rate of key monitoring points in the simulation, the borehole layout correction parameters used to correct the risk of under-excavation are evaluated. The scheme optimization module is used to apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the optimization of the blasting scheme and obtain the optimized tunnel blasting scheme.

[0006] The beneficial effects of this invention are as follows: On the one hand, by directly acquiring the 3D design model of the tunnel with a non-circular contour and identifying the high-risk geometric regions of the tunnel formed by abrupt changes in curvature and normal, a spatial blasting constraint field that can quantify the local geometric complexity is constructed. This process abandons the simplification of the contour by traditional methods and can fully preserve and utilize the crucial spatial geometric information on the contour. Furthermore, by planning the ideal blasting drilling path, it can ensure that the arrangement of blast holes responds preferentially to the constraints of positions with drastic curvature changes, such as the arch waist and arch foot, and achieves adaptive matching between blast hole density and contour complexity. This avoids the blindness of traditional equal-spacing hole layout and can effectively cope with the refined requirements of blasting shaping for complex curved contours, significantly improving the accuracy of blasting scheme design.

[0007] On the other hand, a blasting energy sequence matching the local constraint strength is derived based on the ideal blasting borehole path, and this sequence is used as the boundary condition to perform rock mass damage dynamics blasting simulation. The differentiated allocation of the energy sequence, based on its strong correlation with local geometric risks, enables "on-demand energy supply" in different sections of complex contours, avoiding energy waste and damage. Simultaneously, by tracking the energy transfer attenuation rate of preset simulation monitoring points during the simulation process, not only can the actual effect of blasting energy in the rock mass be dynamically evaluated, but also potential under-excavation risk areas can be accurately predicted and located, significantly improving the accuracy of blasting effect prediction. Attached Figure Description

[0008] Figure 1 This is a flowchart illustrating the steps of the three-dimensional dynamic simulation method for the entire tunnel blasting process of the present invention. Figure 2 This is a block diagram of the three-dimensional dynamic simulation system for the entire tunnel blasting process of the present invention; Figure 3 A three-dimensional digital model of the tunnel structure; Figure 4 A comparative diagram showing the drilling layout optimization before and after; The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0009] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0010] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0011] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0012] To achieve the above objectives, please refer to Figures 1 to 4 This invention provides a three-dimensional dynamic simulation method for the entire process of tunnel blasting, comprising the following steps: Step S1: Obtain a 3D design digital model of the tunnel with a non-circular contour; identify high-risk geometric areas of the tunnel in the 3D design digital model; In this embodiment of the invention, firstly, a 3D tunnel design model file (e.g., DWG or SAT format) designed by AutoCAD Civil 3D or similar software is imported via a data interface. This model precisely defines the tunnel cross-sectional profile, which is composed of a three-centered circle, curved sidewalls, and a straight wall arch. Multiple two-dimensional profile sections are automatically extracted along the tunnel's longitudinal axis at 200 mm intervals, and 3D spatial coordinate points are collected at 100 mm intervals along the profile arc for each section, generating profile spatial point cloud data covering the entire tunnel excavation boundary. Subsequently, based on this point cloud data, target geometric structure regions such as the arch crown, arch waist, and sidewalls are identified by analyzing the point's normal vector and local curvature. Finally, within these target geometric structure regions, three types of high-risk tunnel geometric regions are accurately identified and divided by setting specific geometric thresholds: regions with abrupt changes in normal vector boundaries, regions with abrupt changes in curvature radius, and structural intersection and congestion regions.

[0013] Step S2: Construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric area of ​​the tunnel; plan the ideal blasting borehole path by analyzing the cumulative weight of the blasting constraint values ​​in the spatial blasting constraint field; determine the key points of simulation monitoring according to the high-risk geometric area of ​​the tunnel. In this embodiment of the invention, different risk contribution weights are assigned to three high-risk areas. Then, based on these weights and point cloud data, a continuous three-dimensional scalar field is constructed using a spatial interpolation algorithm. Areas with higher values ​​in this field represent areas with stronger blasting constraints and greater difficulty. Next, based on this constraint field, a planning algorithm similar to "finding the highest cost path" is used to calculate an ideal blasting drilling path that prioritizes traversing all high-constraint areas. Simultaneously, at the geometric center of the high-risk areas, simulated monitoring points for subsequent simulation tracking are automatically deployed.

[0014] Step S3: Obtain the three-dimensional geological model of the tunnel; derive the blasting energy sequence based on the ideal blasting borehole path, and use this as the boundary condition to perform rock mass damage dynamics blasting simulation on the three-dimensional geological model of the tunnel. By tracking the energy transfer attenuation rate of key monitoring points in the simulation, evaluate the borehole layout correction parameters used to correct the risk of under-excavation. In one implementation of this invention, all nodes along the ideal blasting borehole path are traversed. At each node, its constraint value in the spatial blasting constraint field is queried. Assuming the constraint value at the arch waist is 0.8, and the constraint value at the straight wall section is 0.1, the base blasting energy (equivalent TNT) set for the borehole at the arch waist will be 8 times that of the straight wall section, thus forming a blasting energy sequence precisely matched to the constraint strength. The blasting energy sequence is loaded as a boundary condition into the simulation model. After the simulation begins, the peak vibration velocity (PPV) of the simulated monitoring point (e.g., located at the arch waist) is recorded as 1.2 m / s. The vibration velocity decay of the monitoring point within 15 milliseconds after detonation is calculated. The preset time window is 15 milliseconds, and the decay threshold is 80%. If the vibration velocity decays to 0.15 m / s after 15 milliseconds, the energy transfer decay rate is (1.2-0.15) / 1.2≈87.5%, which is greater than the 80% threshold, indicating a risk of under-excavation in the area. Subsequently, the excess portion of the attenuation rate (87.5% - 80% = 7.5%) is converted into a borehole layout correction command according to a preset ratio. The final generated borehole layout correction parameters are: reduce the borehole spacing in the corresponding area by 10% and increase the borehole outboard angle by 1.5 degrees.

[0015] Step S4: Apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the optimization of the blasting scheme and obtain the optimized tunnel blasting scheme.

[0016] In one implementation of this invention, borehole layout correction parameters ("spacing reduced by 10%, external angle increased by 1.5 degrees") are read, and the system automatically locates the section on the ideal blasting borehole path corresponding to the under-excavation risk zone. The spacing of all boreholes on this section is multiplied by 0.9, and the axis vector is rotated 1.5 degrees around the tangent. Furthermore, the energy value at the corresponding position in the blasting energy sequence is also increased by 10%. After the parameter update is completed, verification is performed. If the verification passes, all the finally determined borehole information is encapsulated into a construction data package (e.g., conforming to IFC or a custom XML format) that can be directly read by the intelligent drilling rig. This data package is the final optimized tunnel blasting scheme, guiding precise on-site construction.

[0017] Preferably, identifying high-risk geometric regions of the tunnel in the 3D design model of the tunnel in step S1 includes: Extract multiple two-dimensional design profile sections along the longitudinal axis of the tunnel's three-dimensional design digital model; Along the arc of each two-dimensional design contour section, three-dimensional spatial coordinate points are collected at preset arc length intervals to generate contour spatial point cloud data; The tunnel's geometric topology is analyzed based on the contour spatial point cloud data to identify the target geometric structure regions of the arch, arch waist, and sidewalls. The high-risk geometric regions of the tunnel are determined based on the target geometric structure region; among them, the high-risk geometric regions of the tunnel include the normal abrupt boundary region, the curvature radius abrupt region, and the structural intersection and congestion region.

[0018] In one implementation of this invention, a 3D tunnel design model conforming to the Industrial Basic Class (IFC) standard, generated by 3D computer-aided design software, is loaded. Assuming a tunnel model with a total length of 100 meters is being processed, its vertical axis is a spatial curve. Starting from the beginning of the curve, the model advances along the curve in steps of 200 millimeters. At each step point, a plane perpendicular to the tangent direction at that point is generated, and the intersection of this plane and the 3D tunnel design model is calculated. This intersection is a 2D design profile section. For this 100-meter tunnel, this step will generate 500 discrete 2D design profile sections.

[0019] In this embodiment of the invention, for each two-dimensional design profile section extracted in the previous step, starting from the highest point of its arch, sampling is performed along the arc path of the profile at a preset arc length interval, and the three-dimensional spatial coordinates of each sampling point are accurately recorded. The set of coordinates of all sampling points on all cross sections together constitutes the contour spatial point cloud data describing the entire tunnel excavation boundary.

[0020] It should be noted that the preset arc length interval is determined based on the minimum geometric feature size of the tunnel profile to ensure that critical geometric details are not lost. Specifically, this interval value is usually set to no more than one-tenth of the minimum designed radius of curvature in the tunnel profile, to ensure that there are enough data points for accurate description even in areas with the most dramatic curvature changes.

[0021] In this embodiment of the invention, geometric analysis is performed on the generated contour spatial point cloud data. By judging the positional attributes and local geometric features of each point, it is automatically classified into different target geometric structure regions.

[0022] Specifically, the identification process includes: Arch crown region identification: For each two-dimensional section, the point with the largest Z-axis coordinate value is identified, and all points within a 2-meter arc length on either side of this point are marked as the arch crown; Side wall region identification: The local curvature of each point is calculated. When the local curvature of a continuous point cloud segment is close to 0 (i.e., approximately a straight line), and the Z-axis component of its normal vector is less than 0.1 (i.e., close to horizontal), this segment of the point cloud is marked as a side wall; Arch waist region identification: All remaining point cloud segments not marked as arch crowns or side walls, especially the curvature variation areas connecting the arch crown and side walls, are automatically identified and marked as the arch waist.

[0023] Preferably, determining the high-risk geometric region of the tunnel based on the target geometric structure region includes: The rate of change of the spatial angle between the normal vectors of adjacent points is calculated based on the spatial point cloud data of the contour in the target geometric structure region. Regions with a spatial angle change rate greater than a preset gradient threshold are extracted as normal abrupt boundary regions. Identify all regions in the target geometric structure area whose radius of curvature is less than a preset radius of curvature threshold, forming abrupt changes in radius of curvature. Based on the target geometric structure area, the area where the angle between the arch waist and the side wall is between 85-95° is identified as the structural intersection and congestion area.

[0024] In this embodiment of the invention, a set of points located in the target geometric structure region (such as the transition zone between the arch crown and the arch waist) is first extracted from the contour spatial point cloud data. For any point in this set... Take the two points that are immediately before and after it. and These three points define a tiny plane in space. By calculating the normal vector of this tiny plane, we can obtain the points. Let nᵢ be the unit normal vector at point n. Traverse the set of points, calculate the unit normal vector for each point in turn, and further calculate the unit normal vector for each adjacent point. and Rate of change of the spatial angle between normal vectors .

[0025] It should be noted that the preset gradient threshold is determined based on the requirements for contour smoothness in the smooth blasting acceptance specifications. The higher the smoothness requirement, the smaller the threshold is set to capture more subtle changes in the contour normal. In a specific application scenario of this invention, to address the requirement that high-speed railway tunnels do not allow obvious angles, the threshold is set to 10 degrees.

[0026] In one implementation of this invention, it is assumed that the analysis reaches the connection point between the three-centered circle of the arch crown and the curve of the arch waist. For the point at that location... Calculate its unit normal vector for For its neighboring points Calculate the unit normal vector for The rate of change of the spatial angle between the two can be calculated using the following formula. : ; in, For point Rate of change of spatial angle at that location and Points and The unit normal vector, " represents the vector dot product operation. Substitute the data to calculate: Degrees. Since the calculated degree of 11.7 degrees is greater than the preset gradient threshold of 10 degrees, the point... The region and its adjacent area are marked as the normal mutation boundary region.

[0027] In this embodiment of the invention, the analysis focuses on the point cloud of the target geometric structure region, specifically the arched waist. For any point within this region... Similarly, take the two points that are immediately before and after it. and These three points can uniquely define a circle in space, and the radius of that circle is the radius of the points. Local radius of curvature at The entire arched region will be traversed, and the local radius of curvature at all points will be calculated.

[0028] In one implementation of this invention, three spatial points , , The coordinates are known, and the local radius of curvature is... It can be calculated using Menger's curvature formula.

[0029] Specifically, the preset radius of curvature threshold is determined by considering the minimum operating turning radius of drilling and blasting machinery (such as drilling rigs) and empirical assessment of the stress concentration effect of blasting. In a specific application scenario of this invention, to avoid excessive rock fragmentation due to overly dense blast holes at small radii of curvature, the threshold is set to 5 meters.

[0030] In this embodiment of the invention, the endpoint of the point cloud in the arched waist region and the starting point of the point cloud in the sidewall region are automatically identified, as these two points are topologically adjacent. The best-fit tangent vectors for the endpoint point cloud segments (e.g., the last 5 points) in the arched waist region are calculated respectively. And the direction vector of the cloud segment at the starting point of the sidewall region (e.g., the first 5 points). Then, calculate the spatial angle between the two vectors. This included angle is the angle at which the arch meets the side wall.

[0031] It is important to note that the defined angle identification range is designed to accurately capture structures with "hard turns" approaching 90 degrees. These structures are typical stress concentration points during blasting, and are highly susceptible to over-excavation or root zone damage due to the inability to smoothly transfer energy. The 85-95 degree range is an empirical value determined after summarizing numerous similar engineering blasting examples.

[0032] In one implementation of this invention, it is assumed that the junction area between the arch waist and the sidewall is identified. The calculated tangent vector at the end of the arch waist... for The direction vector of the sidewall for The intersection angle is calculated using the vector dot product formula. : ; Substituting the data, a 45-degree angle is calculated. The result shows this is a gentle transition and not part of a congested area. At another location, assuming the calculated... for , for The calculated intersection angle Adding the 90-degree base angle, the actual intersection angle is approximately 98 degrees, which is also outside the 85-95 degree range.

[0033] Preferably, the construction of a spatial blasting constraint field covering the blasted rock mass based on the high-risk geometric region of the tunnel in step S2 includes: Each three-dimensional spatial coordinate point in the contour spatial point cloud data is defined as a network node, and the node corresponding to the center point of the arch is identified as the starting node. Starting from the initial node, adjacent nodes are connected sequentially along the contour arc based on the topological adjacency relationship of three-dimensional spatial coordinate points to generate a virtual blasting conduction network for the blasting energy transfer path; Based on the high-risk geometric areas of the tunnel, risk contribution weights for the blasting sensitivity of each point in the contour spatial point cloud data are set separately. Based on the risk contribution weight analysis, the blasting constraint values ​​in the virtual blasting transmission network are analyzed to construct a spatial blasting constraint field.

[0034] In this embodiment of the invention, each three-dimensional spatial coordinate point in the contour spatial point cloud data is abstracted as a network node in graph theory. Simultaneously, for the needs of subsequent path analysis, a unique starting point must be determined. By traversing all network nodes, the node with the largest average Z-axis coordinate value among all two-dimensional cross-sections is identified and ultimately determined as the starting node of the entire network. This node physically corresponds to a point on the centerline of the tunnel's arch.

[0035] In this embodiment of the invention, the inherent order relationship (i.e., topological adjacency relationship) of the contour spatial point cloud data during generation is utilized to connect each node to its next node in the data list with a directed edge. This connection process unfolds along the arc path of the contour, from the center of the arch towards one side of the arch foot, until all nodes on one side have been traversed. Returning to the starting node, the same connection operation is performed towards the other side of the arch foot. Ultimately, all nodes and connecting edges together form a chain-like virtual blasting conduction network capable of completely describing the single-loop contour of the tunnel.

[0036] Specifically, the construction of this network does not involve complex geometric calculations, but directly utilizes the spatial order of point cloud data acquisition. For example, points and points They are adjacent in physical space and also adjacent in data storage; they only need to be... and Simply establish a connection between them.

[0037] In this embodiment of the invention, a quantitative risk contribution weight is set to characterize the degree of contribution to the difficulty of blasting formation. This weight is set based on prior knowledge of engineering experience and numerical simulation. For example, the weights are set as follows: Regions with abrupt changes in radius of curvature: Since they directly affect the density and stress concentration of boreholes, they have the greatest impact on formation, and their risk contribution weight is set to 1.5; Regions with abrupt changes in normal direction boundary: These regions are prone to forming sharp angles, affecting the smoothness of the contour, and their risk contribution weight is set to 1.2; Congested areas where structures intersect: These regions are intersections of multiple geometric features and have a high risk, and their risk contribution weight is set to 1.3. For ordinary points that do not belong to any high-risk areas, their weight value is defaulted to 1.0. All nodes in the contour space point cloud data are traversed, and corresponding risk contribution weights are assigned to them according to their respective regions.

[0038] In one implementation of this invention, a three-dimensional mesh is first defined around the tunnel outline, covering the expected blasting influence zone (e.g., 5 meters outside the outline). Then, the blasting constraint value for each three-dimensional mesh point is calculated using inverse distance weighting (IDW) interpolation. The calculation formula is as follows: ; in, The three-dimensional mesh points to be calculated The blasting constraint value at the location; It is the first Risk contribution weight of each contour space point cloud data node (i.e., discrete blasting constraint value). It is the first Each contour node corresponds to a grid point. The weight of is calculated using the following formula: , Grid points With contour nodes The Euclidean distance between them It is a power parameter (usually taken as 2). By performing the above calculations on all three-dimensional mesh points, a continuous three-dimensional scalar field is finally generated that can quantify the difficulty of blasting at any point around the tunnel outline. This field is the spatial blasting constraint field. In this field, the field value is the highest near the region of abrupt change in the radius of curvature, while the field value gradually decreases in the region far from the outline.

[0039] Preferably, the blasting constraint values ​​in the virtual blasting transmission network analyzed according to risk contribution weights include: Obtain the tunnel excavation centerline based on the two-dimensional design profile section; Calculate the shortest spatial distance between each network node in the virtual blasting transmission network and the tunnel excavation centerline; Take the reciprocal of the shortest spatial distance and determine it as the baseline blasting constraint value corresponding to the network node; Determine whether each network node in the virtual blasting transmission network is located within any high-risk geometric region of the tunnel. If the network node is located within a high-risk geometric region of the tunnel, multiply the risk contribution weight corresponding to that region by the node's baseline blasting constraint value, and use the product as the node's blasting constraint value. If a network node is not located in any high-risk geometric area of ​​a tunnel, its baseline blasting constraint value is the blasting constraint value of that node.

[0040] In this embodiment of the invention, each independent two-dimensional cross-section is considered as a closed polygon composed of multiple curves or straight lines. By calculating the geometric centroid of the polygonal region, a center point of the cross-section in space can be obtained. This process is repeated for all two-dimensional design profile cross-sections, calculating the center point of each cross-section. Finally, all these center points are connected sequentially along the tunnel's longitudinal axis to form a three-dimensional spatial curve, which is defined as the tunnel excavation centerline.

[0041] In this embodiment of the invention, each network node in the virtual blasting propagation network will be traversed. For any given network node... This involves calculating the shortest distance from the point to the spatial curve representing the tunnel excavation centerline generated in the previous step. Specifically, this calculation process is an optimization problem of finding the shortest distance from a point to the spatial curve. A point is found on the tunnel excavation centerline. , making the point With point Euclidean distance between Minimum. This minimum distance value. That is, network nodes The shortest spatial distance between the tunnel excavation centerline and the tunnel centerline.

[0042] In this embodiment of the invention, the shortest spatial distance for each network node calculated in the previous step is... A mathematical transformation is performed to generate an initial constraint value. The physical meaning of this transformation is that the farther a point is from the tunnel centerline (usually a sidewall or arch foot), the greater its contribution to the geometric constraints of the profile shaping.

[0043] In this embodiment of the invention, each network node in the virtual blasting conduction network is traversed, and its coordinates are spatially matched with the previously determined high-risk geometric areas of the tunnel. If a network node... If the coordinates fall within a high-risk region (such as a region of abrupt change in curvature radius), the preset risk contribution weight for that region is read. Then, the constraint value of the node is enhanced through multiplication.

[0044] In this embodiment of the invention, for network nodes that, after spatial location matching, do not fall within any high-risk geometric area of ​​the tunnel, they are considered to be constrained only by their geometric location (distance from the centerline) and not affected by local geometric abrupt changes. Therefore, for such nodes, their baseline blasting constraint value is directly used as their final blasting constraint value.

[0045] Preferably, step S2, which involves analyzing the cumulative weighted sum of blasting constraint values ​​in the spatial blasting constraint field to plan the ideal blasting borehole path, includes: Calculate the cumulative weighted sum of path blasting constraint values ​​from the starting node to each node in the network along the topology of the virtual blasting propagation network; Identify the path with the largest cumulative weight and determine that path as the key blasting propagation path.

[0046] In this embodiment of the invention, the virtual blasting conduction network and the final blasting constraint value calculated for each node in the network in the previous step are retrieved. The virtual blasting conduction network is essentially a chain-like graph structure expanding outwards from the starting node (the center of the dome). Starting from the starting node, the network is traversed along its topology (i.e., the node order). During the traversal, the values ​​from the starting node to any currently traversed node are dynamically calculated. The sum of the brute-force constraint values ​​of all nodes along the path. Specifically, this sum is called the cumulative weighted sum, and its calculation formula is as follows: ; in, From the starting node To the current node The cumulative weights sum; It is the first on the path The final blasting constraint value for each node; It is the first on the path Each node and its preceding node The arc length between (i.e., the length of an edge in the virtual blast propagation network). Introducing arc length. As a weighting factor, it is intended to more accurately reflect the cumulative effect of constraints on long-distance paths.

[0047] It should be noted that the virtual blasting propagation network has two branches, corresponding to the left and right halves of the profile, respectively. Therefore, the calculation process is performed independently on the two branches, ultimately resulting in two cumulative weight and variation curves extending from the arch crown to the arch foot.

[0048] In one implementation, assume the starting node Explosion constraint value It is 0.125, which is related to the next node. arc length It is 0.1 meters. Explosion constraint value It is 0.128. Then it reaches... Cumulative weights and The calculation is as follows: (This is a simplified trapezoidal rule approximation for integral calculation), and the calculation will proceed point by point to the endpoint of the contour in this manner.

[0049] In this embodiment of the invention, after calculating the cumulative weights of all nodes in the virtual blasting propagation network, the cumulative weights and total values ​​of the left and right branches of the profile at the final endpoint of the arch foot are compared. These two total values ​​represent the total constraint accumulated from the center of the arch along the left and right profiles to the arch foot. The branch path with the larger cumulative weight and total value is selected as the path with the strongest global constraint and the greatest difficulty in blasting formation.

[0050] Preferably, planning the ideal blasting borehole path by analyzing the cumulative weighted sum of blasting constraint values ​​in the spatial blasting constraint field further includes: Nodes whose blasting constraint values ​​in the spatial blasting constraint field are within a preset screening ratio threshold are selected as core risk nodes. The core risk nodes and key blasting transmission paths are merged and deduplicated to obtain the mandatory path points; Using each point in the forced path as a vertex of the graph, calculate the shortest arc length between any two vertices along the contour spatial point cloud data as the weight of the edge, and construct a weighted risk topology graph. In the weighted risk topology, the path with the minimum cumulative weight that starts from the starting point of the critical blasting transmission path, passes through all forced transit points, and finally reaches the end point of the critical blasting transmission path is defined as the ideal blasting borehole path.

[0051] In this embodiment of the invention, the final blasting constraint values ​​of all network nodes in the spatial blasting constraint field are sorted in descending order. Then, according to a preset screening ratio threshold, a subset of nodes are selected from the top of this sorted list. These selected nodes represent the points with the strongest blasting forming constraints on the entire contour and are collectively marked as a set, namely, the core risk nodes.

[0052] It should be noted that the preset screening ratio threshold is designed to balance computational efficiency and path planning accuracy. This threshold is typically set between 15% and 25% based on engineering experience. Setting it too high will result in too many nodes needing to be processed, increasing the complexity of subsequent path calculations; setting it too low will miss some minor but still important risk points.

[0053] In this embodiment of the invention, the core risk node set extracted in the previous step is combined with the previously determined key blasting transmission path (which is itself a sequence of nodes) using a union operation. The purpose of this operation is to integrate the path with the strongest global constraints with the discrete points with the strongest local constraints. After the operation, a deduplication operation is performed to ensure that each node appears only once in the final set. This new set of nodes obtained after merging and deduplication is defined as the mandatory path point. Specifically, this step ensures that the final planned path not only proceeds along the general direction of the key blasting transmission path with the highest overall difficulty, but also precisely passes through all the core risk nodes with the highest local risk, even if some core risk nodes were not originally on the key blasting transmission path.

[0054] In this embodiment of the invention, each point in the set of forced path points generated in the previous step is abstracted as a vertex in graph theory. Then, the shortest arc length connecting any two vertices along the topological order of the contour space point cloud data is calculated. This shortest arc length is used as the weight of the edge connecting the two vertices. By calculating the shortest arc length between all vertex pairs, a fully connected, undirected, weighted risk topology graph is constructed.

[0055] In this embodiment of the invention, after the weighted risk topology graph is constructed, the goal of the problem is to find a path starting from the starting point of the critical blasting propagation path (which is also one of the mandatory path points), such that the path must visit every other vertex in the set of mandatory path points at least once, and finally reach the ending point of the critical blasting propagation path (which is also one of the mandatory path points), while requiring the total length of the entire path (i.e. the sum of the weights of all the edges traversed) to be minimized.

[0056] Specifically, since this is a computationally complex problem, approximate algorithms, such as genetic algorithms or simulated annealing, are used to solve for the minimum cumulative weight path. The solved path is an ordered sequence containing all forced waypoints. The optimal path sequence obtained after calculation is... Following this sequence, the original paths of these nodes in the contour space point cloud data are connected sequentially to form a continuous, smooth spatial curve. This curve not only follows the general direction with the highest global risk but also precisely connects all the core points with the highest local risk, while ensuring the shortest total path length, thus optimizing the blasting sequence. This curve is ultimately defined as the ideal blasting drilling path.

[0057] Preferably, step S3, deriving the blasting energy sequence based on the ideal blasting borehole path, includes: Based on the spatial blasting constraint field and the virtual blasting conduction network, a three-dimensional minimum resistance line vector field is constructed to characterize the directional difficulty of blasting energy transfer. Select the endpoint in the ideal blasting borehole path and calculate the basic blasting energy required for the terminal blast hole by combining the corresponding vector modulus in the three-dimensional minimum resistance line vector field. The basic blasting energy is calculated iteratively along the ideal blasting borehole path towards the starting point, and a discrete energy point sequence corresponding to the path position is obtained, forming the blasting energy sequence.

[0058] In this embodiment of the invention, a spatial blasting constraint field and a virtual blasting conduction network are loaded. To simulate the path and difficulty of blasting energy transmission from the blast hole to the excavation free face (i.e., the tunnel face), a three-dimensional minimum resistance line vector field is constructed. Specifically, at each node of the virtual blasting conduction network, a three-dimensional vector pointing towards the tunnel excavation free face is defined. The direction of this vector is set to originate from that node, be perpendicular to the tunnel's longitudinal axis, and point towards the interior of the tunnel. The magnitude (i.e., modulus) of the vector is determined by the blasting constraint value of that node in the spatial blasting constraint field.

[0059] In one implementation of this invention, for a node in a virtual blasting conduction network... Its blasting constraint value Assuming the tunnel's longitudinal axis is the Y-axis, then the minimum resistance line vector at this point... It can be represented as: ; in, It is a node The minimum resistance line vector at that point, This is the blast constraint value at that point. It is a unit direction vector perpendicular to the tunnel's longitudinal axis and pointing towards the tunnel's center. By calculating this vector for all nodes in the network, a three-dimensional minimum resistance line vector field covering the entire tunnel outline is finally generated.

[0060] In this embodiment of the invention, the last node, i.e., the endpoint (usually located at the arch foot or the arch crown on the other side), is extracted from the ideal blasting borehole path generated in the previous step. This point represents the last contour point in the entire blasting sequence that requires precise control. The vector corresponding to this endpoint in the three-dimensional minimum resistance line vector field is queried, and its modulus is extracted.

[0061] Specifically, the calculation of basic blast energy is based on classical blast theory formulas, but the key parameter—the minimum resistance line—is replaced by the vector modulus. The calculation formula is as follows: ; in, It is the basic blasting energy required for the final blast hole (unit: kilojoules). It is a coefficient related to rock properties, determined by the rock blastability level in the input three-dimensional geological model; It is the unit consumption coefficient of the explosive, which is determined by the type of explosive selected; It is the modulus of the vector corresponding to the endpoint in the three-dimensional minimum resistance line vector field.

[0062] In this embodiment of the invention, after calculating the basic blasting energy at the endpoint, iterative calculations are performed along the ideal blasting borehole path, starting from the endpoint and proceeding step by step towards the starting point. In each iteration, the energy transfer relationship between the current node and the next node (the node closer to the endpoint) is considered.

[0063] In one implementation of this invention, the iterative calculation logic is as follows: the blasting energy required by the current node is equal to the energy already calculated for the next node, plus the energy attenuation caused by rock resistance between the two. This energy attenuation is proportional to the average value of the three-dimensional minimum resistance vector field modulus between the two nodes. The iterative formula can be expressed as: ; in, The current node The blast energy that needs to be calculated; It is the next node on its path. The calculated blast energy; It is an energy decay coefficient, which is also determined by the properties of the rock; and These are the minimum resistance vector moduli corresponding to the two nodes; It is the arc length between the two nodes. (This is achieved by passing through the endpoint...) Initially, by iterating forward using this formula, the precise blasting energy value required for each node on the ideal blasting borehole path can be calculated. All these nodes and their corresponding energy values ​​together constitute a discrete sequence of energy points, which is the final blasting energy sequence.

[0064] Of particular importance is the construction of a three-dimensional minimum resistance line vector field, based on the spatial blasting constraint field and the virtual blasting conduction network, to characterize the directional difficulty of blasting energy transfer. This includes: The data of the spatial blasting constraint field are arranged logically according to the virtual blasting propagation network to construct discrete constraint sorting data; The initial three-dimensional geological model of the tunnel is spatially meshed to obtain three-dimensional rock mass elements; Based on discrete constraint sorting data, spatial distance inverse interpolation is performed on each three-dimensional rock mass unit to construct a three-dimensional minimum resistance line vector field that characterizes the directional difficulty of blasting energy transfer.

[0065] In this embodiment of the invention, a spatial blasting constraint field and a virtual blasting conduction network are loaded. The spatial blasting constraint field is a dataset in which each node in the virtual blasting conduction network is assigned a blasting constraint value. To facilitate subsequent interpolation calculations, this data is reorganized according to the topological order of the virtual blasting conduction network. Specifically, starting from the initial node of the virtual blasting conduction network, the three-dimensional coordinates of each node and its corresponding blasting constraint value are extracted sequentially along its unique path. This extraction process generates an ordered list, where each item contains the three-dimensional coordinates of a point and a scalar value (constraint value). This ordered list, strictly arranged according to path logic, is defined as discrete constraint sorted data.

[0066] In this embodiment of the invention, an initial three-dimensional geological model of the tunnel is loaded. This model includes not only the cavity geometry after tunnel excavation but also the surrounding rock mass within a certain range. To enable continuous physical field calculations within the rock mass, this continuous rock mass model needs to be discretized. Specifically, finite element preprocessing technology is used, employing tetrahedral or hexahedral elements to spatially mesh the entire rock mass model. After meshing, the original continuous rock mass model is transformed into a collection of numerous closely connected, minute three-dimensional rock mass elements. Each three-dimensional rock mass element possesses its own volume, center coordinates, and rock mass mechanical properties inherited from the geological model.

[0067] In this embodiment of the invention, all three-dimensional rock mass elements generated in the previous step are traversed. For each three-dimensional rock mass element, its geometric center point is taken as the representative point of that element. Then, based on the discrete constraint sorting data, an interpolated blasting constraint value is calculated for this representative point using the inverse distance weighting (IDW) method. This interpolation process actually smoothly extends the constraint information discretely distributed on the contour line to the entire three-dimensional rock mass space.

[0068] Preferably, in step S3, the blasting energy sequence is derived based on the ideal blasting borehole path, and this sequence is used as a boundary condition to perform rock mass damage dynamics blasting simulation on the three-dimensional geological model of the tunnel. By tracking the energy transfer attenuation rate of key monitoring points in the simulation, the borehole layout correction parameters used to correct the risk of under-excavation are evaluated, including: The blasting energy sequence is used as the energy boundary condition, and the ideal blasting borehole path is used as the detonation sequence guide to construct the blasting simulation task. Based on the blasting simulation task, rock mass damage dynamics blasting simulation was performed on the three-dimensional geological model of the tunnel to obtain the simulated blasting process; During the simulated blasting process, data is tracked at all key locations marked by the simulation monitoring to obtain the critical simulation timeline. The peak vibration velocity of each simulated monitoring point in the key simulation time history is analyzed, and its attenuation rate within the preset time window after detonation is calculated to obtain the energy transfer attenuation rate. When the energy transfer attenuation rate is greater than the preset attenuation threshold, the simulated monitoring focus is marked as an area with under-excavation risk, and the process is traced back along the ideal blasting borehole path to determine the main blasting conduction path that provides energy. The average spacing of blast holes along the main blasting propagation path is calculated and compared with the average spacing of blast holes in areas with under-excavation risk. The hole spacing adjustment coefficient and the drilling angle adjustment amount are dynamically evaluated and together constitute the drilling layout correction parameters.

[0069] In this embodiment of the invention, previously generated core data—the blasting energy sequence and the ideal blasting borehole path—are integrated. Specifically, the ideal blasting borehole path defines the precise spatial location of the boreholes and the order of detonation, while the blasting energy sequence assigns a specific energy value (equivalent to TNT) to each borehole location along this path. This data is then converted into a script file conforming to the input format of dynamic simulation software (such as LS-DYNA). In this script file, each borehole location is defined as a blasting load application point, its energy-time curve determined by the energy value in the blasting energy sequence, and the detonation delay time determined by its order within the ideal blasting borehole path. This script file, containing all load definitions, loading order, and simulation control parameters, constitutes a complete blasting simulation task.

[0070] In this embodiment of the invention, the blasting simulation task constructed in the previous step is submitted to a dynamics solver deployed on a high-performance computing server. The solver loads a three-dimensional geological model of the tunnel (a finite element model containing millions of rock mass elements) and applies the boundary conditions and loads defined in the blasting simulation task. It should be noted that the simulation uses a constitutive model capable of describing the damage and fracture behavior of rock under high strain rates, such as the HJC (Holmquist-Johnson-Cook) model. The solver performs explicit dynamic calculations in microsecond-level time steps, simulating the entire process from the moment of detonation, including the propagation and reflection of the shock wave in the rock mass, and the resulting accumulation of damage to rock mass elements and the propagation of macroscopic cracks. This complete sequence of calculation results, lasting from several milliseconds to tens of milliseconds, constitutes the simulated blasting process.

[0071] In this embodiment of the invention, data probes previously set up at key monitoring points (e.g., high-risk locations at the arch waist and arch foot) are activated simultaneously with the start of the simulation calculation. Throughout the simulated blasting process, these data probes record the changes in various physical quantities of the three-dimensional rock mass elements at their locations in real time with extremely high temporal resolution (e.g., every 0.1 microseconds), particularly velocity, acceleration, and stress state. All these time-varying data are continuously recorded, forming a series of time-series data files, which constitute the key simulation time history.

[0072] In one implementation of this invention, the preset time window is set according to the brittle characteristics of the rock and the blasting design requirements, for example, 15 milliseconds. Assuming the simulated monitoring point located at the arch waist has a peak vibration velocity of 1.2 m / s, occurring 5.2 milliseconds after detonation, the vibration velocity at 20.2 milliseconds (5.2 + 15) after detonation will be recorded, assumed to be 0.15 m / s. Then the energy transfer attenuation rate ( The calculation formula for ) is as follows: ; in, It is the peak vibration velocity. It is the speed at the end of the time window; substitute the data to calculate: ;this The energy transfer attenuation rate is the assessment result of this monitoring focus.

[0073] In this embodiment of the invention, the energy transfer attenuation rate calculated in the previous step is compared with a preset attenuation threshold. This threshold is a critical value determined based on engineering experience to determine whether the rock can be effectively broken, and is usually set between 75% and 85%. If the energy transfer attenuation rate exceeds this threshold, it means that the energy attenuates too quickly in this area and fails to cause sufficient and sustained damage to the rock. Therefore, the area where the simulated monitoring focus is located is marked as an area with under-excavation risk. Subsequently, along the ideal blasting borehole path, starting from the location of the monitoring point, a reverse search is performed towards the starting point of the detonation to find the sequence of blast holes closest to the monitoring point. This sequence is then determined as the main blasting conduction path that plays a major role in breaking the rock in this area.

[0074] In another implementation of this invention, for the under-excavation risk area marked in the previous step: The hole spacing adjustment coefficient is evaluated: the excess portion of the energy transfer attenuation rate in this area (87.5% - 80% = 7.5%) is converted into a hole spacing reduction ratio proportional to the excess through a preset nonlinear mapping function. For example, a 7.5% excess attenuation rate is mapped to a 10% spacing reduction. This 10% is defined as the hole spacing adjustment coefficient. The drilling angle adjustment amount is evaluated: simultaneously, in the final contour shape of the simulation results, the point furthest from the design contour within the under-excavation risk area (i.e., the "root point") is identified. Then, the direction vector pointing from the borehole axis closest to the root point on the main blasting transmission path to the root point is calculated. The spatial angle between this direction vector and the current borehole axis vector is set as the drilling angle adjustment amount used to correct the direction of blasting energy. The hole spacing adjustment factor (a 10% reduction) and the drilling angle adjustment amount (a 1.5-degree increase) are integrated to form a specific and executable set of drilling layout correction parameters.

[0075] Of particular importance, the execution of the simulated blasting process includes: The energy value in the blasting energy sequence is applied to the blasting hole element corresponding to the ideal blasting drilling path as the initial load; The expansion process of detonation products in the detonation borehole was simulated using the smooth particle hydrodynamics method, and the calculated borehole wall pressure was used as the dynamic boundary condition. The evolution of damage variables in each three-dimensional rock mass element in a three-dimensional geological model of a tunnel under dynamic boundary conditions is calculated based on the theory of damage mechanics of continuous media. When the damage variable of a three-dimensional rock mass element reaches 1, the element is determined to be completely destroyed, forming cracks. The calculation continues iteratively until all blasting energy is dissipated, and the simulated blasting process is finally obtained.

[0076] In this embodiment of the invention, within the finite element mesh of the tunnel's three-dimensional geological model, a set of three-dimensional rock mass elements corresponding to the coordinates of each node on the ideal blasting borehole path is precisely located. These elements collectively constitute a virtual blast hole. Then, the blasting energy sequence is read. Specifically, for the first node on the ideal blasting borehole path... Each node corresponds to an energy value in the blast energy sequence. Through an internal energy loading function, it is uniformly applied to the components constituting the first... This process is applied to all three-dimensional rock mass elements of each blast hole. It is completed at time step 0 of the simulation, instantaneously converting chemical energy into the initial internal energy of the material within the blast hole, serving as the initial load for the entire dynamic simulation.

[0077] In this embodiment of the invention, after applying the initial load, a special numerical method suitable for simulating the expansion of high-temperature and high-pressure gases—Smoothed Particle Hydrodynamics (SPH)—is used to accurately simulate the complex behavior of detonation products within the detonation borehole. It should be noted that the SPH method discretizes the detonation products into a series of particles carrying physical properties (such as mass, density, and pressure). These particles move and interact according to the Jones-Wilkins-Lee (JWL) equation of state, which precisely describes the pressure change with volume after the explosive detonation.

[0078] Specifically, the SPH solver calculates the pressure exerted by each particle on surrounding particles at each time step. When these particles move and collide with the borehole wall (i.e., the interface with ordinary 3D rock mass elements), the instantaneous pressure applied to the borehole wall is precisely calculated. This rapidly changing borehole wall pressure is used as the dynamic boundary condition driving rock mass failure and is applied in real time to the 3D rock mass elements in contact with it.

[0079] In this embodiment of the invention, after obtaining the dynamic boundary condition of borehole wall pressure, the mechanical response of all three-dimensional rock mass elements in the entire three-dimensional geological model of the tunnel is calculated. Specifically, a theoretical framework based on continuum damage mechanics is adopted. Under this framework, each three-dimensional rock mass element is introduced with an internal state variable called the damage variable D. The value of this variable ranges from 0 to 1, where D=0 indicates that the material is intact and D=1 indicates that the material has completely lost its load-bearing capacity.

[0080] In this embodiment of the invention, the damage variable D of each three-dimensional rock mass element is monitored in real time during continuous iterative calculations. When the damage variable D of a certain element accumulates and reaches the critical value of 1 for the first time, the element is immediately determined to be completely destroyed and has lost all mechanical strength. In the visualization post-processing, this element judged to be completely destroyed can be rendered as part of a crack or directly deleted from the model, thereby intuitively showing the formation and propagation path of macroscopic cracks.

[0081] Of particular importance is that, based on the high-risk geometric areas of the tunnel, the focus of simulated monitoring includes: K-means clustering was performed on the risk values ​​of high-risk geometric areas of the tunnel to identify areas prone to over- or under-excavation during blasting, thus obtaining the core risk area cluster. Identify the geometric center of each cluster within the core risk area cluster to determine the initial simulation monitoring points; Based on the gradient direction of risk values ​​in the spatial blasting constraint field, the sensitivity of the core risk region cluster on the energy transfer path is evaluated, and then the key monitoring points in the simulation process are selected and determined.

[0082] In this embodiment of the invention, all contour spatial point cloud data nodes within the high-risk geometric region of the tunnel, as well as the corresponding blasting constraint values ​​of these nodes in the spatial blasting constraint field, are extracted. These constraint values ​​are considered as risk quantification indicators for each point. To identify areas where risk is highly concentrated in space, an unsupervised machine learning algorithm—K-means clustering—is employed.

[0083] Specifically, the K-means algorithm uses the coordinates of these high-risk points in three-dimensional space as input features and iteratively groups them. The goal of the algorithm is to find K cluster centers that minimize the sum of the squared distances from all points to their respective cluster centers. It should be noted that the number of clusters, K, is preset based on the number of typical risk areas in the tunnel cross-section. For example, for a typical three-center circular cross-section, its main risk points are usually concentrated in 3 to 4 areas, such as the left arch waist, right arch waist, and arch foot; therefore, the K value can be preset to 4. After clustering, the set of point clouds in each category constitutes a core risk area cluster that is spatially adjacent and has similar risk values.

[0084] In this embodiment of the invention, after obtaining each core risk area cluster, a most representative monitoring location is determined for each cluster. Specifically, the arithmetic mean of the three-dimensional coordinates of all point cloud nodes within each core risk area cluster is calculated. This average point is the geometric center of the cluster. The geometric center points of all core risk area clusters are collectively determined as a set, namely the initial simulation monitoring points. These points physically represent the central locations of various high-incidence areas of blasting over-excavation and under-excavation.

[0085] In this embodiment of the invention, not all initial simulation monitoring points are of equal importance. To select the key points that are most sensitive to changes in blast energy and best reflect the effectiveness of the blast, the sensitivity assessment is based on gradient analysis of the spatial blast constraint field. The gradient is a vector pointing in the direction of the fastest growth of the scalar field (in this case, the risk value), and its magnitude represents the rate of growth.

[0086] The present invention also provides a three-dimensional dynamic simulation system 100 for the entire process of tunnel blasting, which executes the three-dimensional dynamic simulation method for the entire process of tunnel blasting as described above. The three-dimensional dynamic simulation system for the entire process of tunnel blasting includes: The geometric risk identification module 101 is used to acquire a three-dimensional design digital model of a tunnel with a non-circular contour; and to identify high-risk geometric areas of the tunnel in the three-dimensional design digital model. The borehole path planning module 102 is used to construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric area of ​​the tunnel; by analyzing the cumulative weight of the blasting constraint values ​​in the spatial blasting constraint field, it plans the ideal blasting borehole path; and determines the simulation monitoring focus according to the high-risk geometric area of ​​the tunnel. The blasting effect simulation module 103 is used to acquire a three-dimensional geological model of the tunnel; blasting energy sequence is derived based on the ideal blasting borehole path, and rock mass damage dynamics blasting simulation is performed on the three-dimensional geological model of the tunnel using this as boundary condition; by tracking the energy transfer attenuation rate of key monitoring points, the borehole layout correction parameters used to correct the risk of under-excavation are evaluated. The scheme optimization module 104 is used to apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the blasting scheme optimization, and obtain the optimized tunnel blasting scheme.

[0087] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0088] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A three-dimensional dynamic simulation method for the entire process of tunnel blasting, characterized in that, Includes the following steps: Step S1: Obtain the 3D design digital model of the tunnel with a non-circular contour; Identify high-risk geometric regions in the 3D design model of a tunnel; Step S2: Construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric region of the tunnel, including: setting the risk contribution weight of the blasting sensitivity of each point in the contour spatial point cloud data according to the high-risk geometric region of the tunnel; Based on the risk contribution weight analysis, the blasting constraint values ​​in the virtual blasting propagation network are analyzed to construct a spatial blasting constraint field, specifically including: Obtain the tunnel excavation centerline based on the two-dimensional design profile section; Calculate the shortest spatial distance between each network node in the virtual blasting transmission network and the tunnel excavation centerline; Take the reciprocal of the shortest spatial distance and determine it as the baseline blasting constraint value corresponding to the network node; Determine whether each network node in the virtual blasting transmission network is located within any high-risk geometric region of the tunnel. If the network node is located within a high-risk geometric region of the tunnel, multiply the risk contribution weight corresponding to that region by the node's baseline blasting constraint value, and use the product as the node's blasting constraint value. If a network node is not located in any high-risk geometric area of ​​a tunnel, its baseline blasting constraint value is the blasting constraint value of that node. By analyzing the cumulative weighted sum of blasting constraint values ​​in the spatial blasting constraint field, an ideal blasting borehole path is planned, including: Calculate the cumulative weighted sum of path blasting constraint values ​​from the starting node to each node in the network along the topology of the virtual blasting propagation network; Identifying the path with the largest cumulative weight sum and determining it as the critical blasting propagation path; planning the ideal blasting borehole path by analyzing the cumulative weight sum of blasting constraint values ​​in the spatial blasting constraint field also includes: Nodes whose blasting constraint values ​​in the spatial blasting constraint field are within a preset screening ratio threshold are selected as core risk nodes. The core risk nodes and key blasting transmission paths are merged and deduplicated to obtain the mandatory path points; Using each point in the forced path as a vertex of the graph, calculate the shortest arc length between any two vertices along the contour spatial point cloud data as the weight of the edge, and construct a weighted risk topology graph. In the weighted risk topology map, the path with the minimum cumulative weight, starting from the starting point of the critical blasting transmission path, passing through all forced transit points, and finally reaching the end point of the critical blasting transmission path, is defined as the ideal blasting borehole path; the focus of simulation monitoring is determined based on the high-risk geometric areas of the tunnel. Step S3: Obtain the three-dimensional geological model of the tunnel; derive the blasting energy sequence based on the ideal blasting borehole path, including: Based on the spatial blasting constraint field and the virtual blasting conduction network, a three-dimensional minimum resistance line vector field is constructed to characterize the directional difficulty of blasting energy transfer. Select the endpoint in the ideal blasting borehole path and calculate the basic blasting energy required for the terminal blast hole by combining the corresponding vector modulus in the three-dimensional minimum resistance line vector field. Based on the basic blasting energy, iterative calculations are performed step by step along the ideal blasting borehole path towards the starting point to obtain a discrete energy point sequence corresponding to the path position, forming a blasting energy sequence. This sequence is then used as a boundary condition to perform rock mass damage dynamics blasting simulation on the tunnel's three-dimensional geological model. By tracking and simulating the energy transfer attenuation rate of key monitoring points, borehole layout correction parameters used to correct under-excavation risks are evaluated. Step S4: Apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the optimization of the blasting scheme and obtain the optimized tunnel blasting scheme.

2. The three-dimensional dynamic simulation method for the entire tunnel blasting process according to claim 1, characterized in that, Step S1 involves identifying high-risk geometric regions in the tunnel's 3D design model, including: Extract multiple two-dimensional design profile sections along the longitudinal axis of the tunnel's three-dimensional design digital model; Along the arc of each two-dimensional design contour section, three-dimensional spatial coordinate points are collected at preset arc length intervals to generate contour spatial point cloud data; The tunnel's geometric topology is analyzed based on the contour spatial point cloud data to identify the target geometric structure regions of the arch, arch waist, and sidewalls. The high-risk geometric regions of the tunnel are determined based on the target geometric structure region; among them, the high-risk geometric regions of the tunnel include the normal abrupt boundary region, the curvature radius abrupt region, and the structural intersection and congestion region.

3. The three-dimensional dynamic simulation method for the entire tunnel blasting process according to claim 2, characterized in that, Based on the target geometric structure region, the high-risk geometric regions of the tunnel include: The rate of change of the spatial angle between the normal vectors of adjacent points is calculated based on the spatial point cloud data of the contour in the target geometric structure region. Regions with a spatial angle change rate greater than a preset gradient threshold are extracted as normal abrupt boundary regions. Identify all regions in the target geometric structure area whose radius of curvature is less than a preset radius of curvature threshold, forming abrupt changes in radius of curvature. Based on the target geometric structure area, the area where the angle between the arch waist and the side wall is between 85-95° is identified as the structural intersection and congestion area.

4. The three-dimensional dynamic simulation method for the entire tunnel blasting process according to claim 2, characterized in that, Step S2 involves constructing a spatial blasting constraint field covering the blasted rock mass based on the high-risk geometric region of the tunnel, including: Each three-dimensional spatial coordinate point in the contour spatial point cloud data is defined as a network node, and the node corresponding to the center point of the arch is identified as the starting node. Starting from the initial node, adjacent nodes are connected sequentially along the contour arc based on the topological adjacency relationship of three-dimensional spatial coordinate points to generate a virtual blasting conduction network for the blasting energy transfer path; Based on the high-risk geometric areas of the tunnel, risk contribution weights for the blasting sensitivity of each point in the contour spatial point cloud data are set separately. Based on the risk contribution weight analysis, the blasting constraint values ​​in the virtual blasting transmission network are analyzed to construct a spatial blasting constraint field.

5. The three-dimensional dynamic simulation method for the entire tunnel blasting process according to claim 1, characterized in that, In step S3, the blasting energy sequence is derived based on the ideal blasting borehole path, and this sequence is used as a boundary condition to perform a rock mass damage dynamics blasting simulation on the three-dimensional geological model of the tunnel. By tracking and monitoring the energy transfer attenuation rate of key areas in the simulation, the borehole layout correction parameters used to correct the risk of under-excavation are evaluated, including: The blasting energy sequence is used as the energy boundary condition, and the ideal blasting borehole path is used as the detonation sequence guide to construct the blasting simulation task. Based on the blasting simulation task, rock mass damage dynamics blasting simulation was performed on the three-dimensional geological model of the tunnel to obtain the simulated blasting process; During the simulated blasting process, data is tracked at all key locations marked by the simulation monitoring to obtain the critical simulation timeline. The peak vibration velocity of each simulated monitoring point in the key simulation time history is analyzed, and its attenuation rate within the preset time window after detonation is calculated to obtain the energy transfer attenuation rate. When the energy transfer attenuation rate is greater than the preset attenuation threshold, the simulated monitoring focus is marked as an area with under-excavation risk, and the process is traced back along the ideal blasting borehole path to determine the main blasting conduction path that provides energy. The average spacing of blast holes along the main blasting propagation path is calculated and compared with the average spacing of blast holes in areas with under-excavation risk. The hole spacing adjustment coefficient and the drilling angle adjustment amount are dynamically evaluated and together constitute the drilling layout correction parameters.

6. A three-dimensional dynamic simulation system for the entire process of tunnel blasting, characterized in that, For executing the three-dimensional dynamic simulation method for the entire tunnel blasting process as described in claim 1, the three-dimensional dynamic simulation system for the entire tunnel blasting process includes: The geometric risk identification module is used to acquire a 3D design model of a tunnel with a non-circular contour; and to identify high-risk geometric areas in the 3D design model of the tunnel. The borehole path planning module is used to construct a spatial blasting constraint field covering the blasting rock mass based on the high-risk geometric area of ​​the tunnel; by analyzing the cumulative weight of the blasting constraint values ​​in the spatial blasting constraint field, it plans the ideal blasting borehole path; and determines the key points of simulation monitoring based on the high-risk geometric area of ​​the tunnel. The blasting effect simulation module is used to acquire a three-dimensional geological model of the tunnel; blasting energy sequence is derived based on the ideal blasting borehole path, and rock mass damage dynamics blasting simulation is performed on the three-dimensional geological model of the tunnel using this as boundary condition; by tracking the energy transfer attenuation rate of key monitoring points in the simulation, the borehole layout correction parameters used to correct the risk of under-excavation are evaluated. The scheme optimization module is used to apply the borehole layout correction parameters to the preset tunnel blasting scheme, update the spacing, angle and charge amount of the blast holes to complete the optimization of the blasting scheme and obtain the optimized tunnel blasting scheme.

Images