A natural caving method rock mass space-time evolution simulation method

Abstract

This invention relates to the field of simulation technology, specifically to a method for simulating the spatiotemporal evolution of rock masses using the natural collapse method. The method includes the following steps: extracting joint surface orientation and intersection points, screening closed path lines, analyzing the direction of deformation rate difference, comparing slip and contact directions to extract dominant nodes, determining the formation of extended boundaries in intersection areas, integrating the time series of slip segments, and generating the spatiotemporal evolution results of the rock mass using the natural collapse method. In this invention, closed path lines are constructed through joint surface orientation and intersection relationships to enhance the complete expression of the structural configuration. The temporal continuity of the deformation rate difference direction is combined to identify continuously changing areas, improving the dynamic tracking capability of deformation trends. Dominant node groups are extracted based on the differences in slip and contact directions. Extended boundaries are reorganized using the path intersection distribution. The time sequence and displacement direction of slip segments are integrated to construct trajectory chains, achieving fine reconstruction of the structural configuration, accurate capture of trend changes, and continuous expression of the evolution trajectory.

Description

Technical Field

[0001] This invention relates to the field of simulation technology, and in particular to a method for simulating the spatiotemporal evolution of rock masses using the natural collapse method. Background Technology

[0002] The field of simulation technology primarily involves the digital modeling and operational prediction of real-world physical systems, natural processes, or engineering structures. It utilizes computer simulation methods to dynamically analyze the behavior of these systems under different environmental, temporal, or operational conditions. Core aspects of this technology include the construction of modeling methods, the setup of numerical computation processes, the handling of the evolution of spatiotemporal variables, and the control of simulation accuracy and stability. It is widely used in the predictive analysis and design verification of complex systems in geological engineering, civil engineering, energy development, and aerospace. Its overall technological development emphasizes improving model accuracy, optimizing computational efficiency, and verifying the consistency between simulation results and actual observations, relying on powerful computing resources and multiphysics coupling modeling capabilities. Among them, the traditional natural collapse method for simulating the spatiotemporal evolution of rock masses refers to a method in rock engineering that uses the natural collapse method to simulate the evolution behavior of rock masses in different time and space ranges. The technical issues it addresses are modeling and reproducing the collapse process and evolution path of rock mass structures under gravity or geological stress release conditions without external force intervention. Traditional methods usually use the discrete element method to divide the rock mass into particles, and based on stress analysis, set contact relationships and failure criteria. They trigger the collapse behavior by using static loading or self-weight reconstruction, and iteratively update the displacement field and velocity field in the calculation. At the same time, the collapse evolution process is controlled by setting the time step to realize the simulation of the rock mass structure state at different time points.

[0003] Traditional natural collapse rock mass evolution simulation uses discrete element method to divide particles and triggers structural response based on stress conditions. However, it does not establish a systematic analytical mechanism for the spatial orientation and intersection relationship of joint surfaces. This results in the inability to form a closed geometric expression of the structural configuration. In deformation behavior modeling, the continuity analysis of temporal difference is not introduced, making it difficult to accurately identify the continuous evolution path of regional deformation trends. The determination of slip node direction does not incorporate the relative relationship between contact directions, and the intersection distribution of path lines is not effectively used for boundary determination of the propagation area. As a result, the structural evolution results are difficult to form a clear and continuous spatiotemporal trajectory expression, affecting the continuity of the overall simulation and the accuracy of the reconstruction of the structural evolution trend. Summary of the Invention

[0004] To achieve the above objectives, the present invention adopts the following technical solution: a method for simulating the spatiotemporal evolution of rock masses using the natural collapse method, comprising the following steps:

[0005] S1: Extract the orientation of the three-dimensional joint surfaces of the rock mass, obtain the location of the intersection points, compare the spatial distance and direction between the intersection points, filter the closable path lines, sort them according to spatial continuity and splice the paths to generate joint closure configuration blocks;

[0006] S2: Based on the joint closure configuration block, extract the regional deformation rate, compare the rate difference and direction of adjacent regions, filter regions with continuous difference direction, and generate deformation trend association blocks;

[0007] S3: Based on the deformation trend association block, extract the shear surface node, compare the difference between the node slip and contact direction, filter the node group with the same direction and sort them into line segments, determine the slip propulsion direction, and generate the slip path main control block;

[0008] S4: Extract the path intersection positions based on the sliding path master control block, count the number of intersections and verify the intersection area range, determine the concentrated propagation area, reorganize to form an extended structure outline, and generate a stable boundary extended structure block;

[0009] S5: Extract slip segments from the extended structure blocks of the stable boundary, arrange them in time sequence and record the displacement direction, and combine the line segment blocks into slip trajectories in time sequence to generate the spatiotemporal evolution results of rock mass by natural collapse method.

[0010] As a further embodiment of the present invention, the joint closure configuration block includes spatial distribution of intersection points, network of closed path lines, and configuration boundary outline; the deformation trend association block includes deformation rate difference zone, direction change trend zone, and continuously evolving boundary circle; the slip path master control block includes consistent slip node group, slip direction dominant line, and node connection sequence line segment; the stable boundary extension structure block includes path intersection central area, propagation channel integration area, and extended configuration boundary; and the spatiotemporal evolution result of the rock mass obtained by the natural collapse method includes time series slip trajectory, directional displacement annotation group, and spatiotemporal composite evolution graph.

[0011] As a further aspect of the present invention, the region where the difference direction is continuous refers to a region where the direction of the difference in deformation rate between adjacent regions is basically consistent and continuously distributed in space.

[0012] As a further aspect of the present invention, the concentrated propagation area refers to the area where sliding paths frequently intersect in space and where sliding activities occur in a highly concentrated manner.

[0013] As a further aspect of the present invention, the specific steps of S1 are as follows:

[0014] S101: Based on the spatial orientation of joint surfaces in the three-dimensional structure of rock mass, extract the intersection coordinates of intersecting joint surfaces, calculate the spatial distance and orientation angle between intersection points, and filter them by combining preset angle and distance thresholds to generate a sequence of spatial positions of intersection points;

[0015] S102: Call the spatial location sequence of the intersection point, calculate the difference between the Euclidean distance and the direction, remove inconsistent point pairs according to the spatial closure threshold and the included angle threshold, construct an adjacency matrix that meets the connection conditions, and generate a node map of joint closable paths.

[0016] S103: Based on the node map of the joint closable path, extract the closed path and sort it according to spatial continuity, connect the node coordinates in the path, splice them to form a closed boundary contour, and generate a joint closure configuration block.

[0017] As a further aspect of the present invention, the specific steps of S2 are as follows:

[0018] S201: Based on the joint closure configuration block, extract the deformation monitoring coordinates of the corresponding region, collect the displacement values ​​of the coordinate points at different times, calculate the deformation rate per unit time, and generate a regional deformation rate matrix.

[0019] S202: Call the region deformation rate matrix, calculate the rate difference between adjacent points and obtain the direction vector, record the points where the direction difference remains consistent in multiple time series, filter the set of continuously shifting points, and obtain the index set of continuously changing regions.

[0020] S203: Based on the continuously changing region index set, establish spatial adjacency relationships and construct closed boundaries, connect coordinates to form contour boundaries, and generate deformation trend association blocks.

[0021] As a further aspect of the present invention, the specific steps of S3 are as follows:

[0022] S301: Based on the deformation trend association block, obtain the node position of the shear surface in the area, compare the angle between the node sliding direction and the contact direction, and filter the node group with the same direction according to the direction difference threshold. Arrange the nodes into a linear sequence according to the connection order to generate a node sequence with the same direction.

[0023] S302: Call the sequence of nodes with the same direction, calculate the direction vector of adjacent nodes, and determine whether they constitute a valid line segment based on the sliding advance threshold. The line segments that meet the advance conditions are grouped into a set according to the node order to obtain the sliding advance line segment set.

[0024] S303: Based on the set of sliding propulsion line segments, retrieve the line segment direction vectors and aggregate them according to the consistency of direction, gather the continuous line segments into a single path structure, and generate the sliding path master control block.

[0025] As a further aspect of the present invention, the specific steps of S4 are as follows:

[0026] S401: Based on the sliding path master control block, extract the intersection coordinates between all path lines in the path line set, count the distribution frequency of the intersection coordinates in the block, and construct a spatial distribution heat map based on the density of coordinates. Extract the contours of continuous clustered areas in the heat map to generate the spatial distribution range of the intersection area.

[0027] S402: Call the spatial distribution range of the intersection area, determine the number of overlapping paths and the density of intersection points within the intersection area, mark the area with an overlap rate exceeding the path overlap threshold as a concentrated propagation area, and extract the boundary in combination with the path extension direction to obtain the boundary of the concentrated propagation area of ​​the path.

[0028] S403: Based on the boundary of the path-centralized propagation area, the coordinates of control points on the contour line of the aggregated area are used to arrange the directions and reconstruct the node connections according to the principle of consistency of direction, forming a closed structural framework and performing boundary fitting to generate stable boundary extension structural blocks.

[0029] As a further aspect of the present invention, the specific steps of S5 are as follows:

[0030] S501: Based on the stable boundary extended structure block, extract the extended structure outline within the block, call the time series labels of the sliding line segments within the outline, rearrange the line segments according to the time progression order, and record the spatial displacement direction corresponding to the line segments to generate a time series sliding line segment set;

[0031] S502: Call the time series sliding line segment set, merge the corresponding blocks of the line segments in the same contour in time order, and maintain the consistency constraint on the displacement direction of each line segment. Connect the sequence blocks to construct a continuous sliding path to obtain the sliding trajectory sequence structure.

[0032] S503: Based on the slip trajectory sequence structure, extract the spatial boundaries and time labels of the trajectory blocks, establish the spatial superposition relationship of the trajectories in different time periods, and generate the spatiotemporal evolution results of the rock mass by the natural collapse method.

[0033] As a further aspect of the present invention, the continuous sliding path refers to a sequence of sliding line segments that are continuous in time sequence, have the same displacement direction in space, and are formed by piecing together blocks.

[0034] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0035] In this invention, closed path lines are constructed by the orientation and intersection of joint surfaces to enhance the complete expression of the structural configuration. By combining the temporal continuity of the deformation rate difference direction, continuously changing regions are identified, improving the dynamic tracking capability of deformation trends. Based on the difference between slip and contact directions, dominant node groups are extracted. The path intersection distribution is used to reorganize and extend the boundary. The temporal sequence and displacement direction of slip segments are integrated to construct a trajectory chain, thereby achieving fine reconstruction of the structural configuration, accurate capture of trend changes, and continuous expression of evolution trajectory. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 This is a schematic diagram of the steps of the present invention;

[0038] Figure 2 This is a detailed schematic diagram of S1 of the present invention;

[0039] Figure 3 This is a detailed schematic diagram of S2 of the present invention;

[0040] Figure 4 This is a detailed schematic diagram of S3 of the present invention;

[0041] Figure 5 This is a detailed schematic diagram of S4 of the present invention;

[0042] Figure 6 This is a detailed schematic diagram of S5 of the present invention. Detailed Implementation

[0043] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0044] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0045] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, their intended meanings are consistent. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, their intended meanings are consistent.

[0046] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0047] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0048] Please see Figure 1 This invention provides a method for simulating the spatiotemporal evolution of rock masses using the natural collapse method, comprising the following steps:

[0049] S1: Based on the spatial orientation of joint surfaces in the three-dimensional structure of the rock mass, the location of the intersection point is obtained. The spatial distance and directional differences between the intersection points are compared. Based on the comparison results, the path lines that can form a closed orientation are selected. The path line position sequence is called, sorted according to spatial continuity, and spliced ​​into a closed boundary contour to generate a joint closure configuration block.

[0050] S2: Based on the joint closure configuration block, extract the deformation rate value of the corresponding location in the region range, perform comparison on the rate difference between adjacent regions and record the direction of the difference, filter the continuously changing region based on the continuity of the difference direction in the time series, call the selected region location to construct the boundary of the continuous deformation area, and generate deformation trend association block;

[0051] S3: Based on the deformation trend association block, the shear surface node position within the area is called through the area, and the direction difference comparison between the sliding direction and the contact direction of the node is performed. Based on the comparison result, the node group with the same direction is selected and arranged into continuous line segments according to the node connection order. The line segments are called to determine the sliding advancement direction and are gathered into the dominant path to generate the sliding path master control block.

[0052] S4: Based on the sliding path master control block, call the path line set to call the intersection position between the path lines, perform statistics on the number distribution of the intersection position and record the intersection area range, determine the concentrated propagation area between the path lines based on the intersection area range, call the propagation area position to reorganize into the extended structure outline, and generate a stable boundary extended structure block;

[0053] S5: Based on the stable boundary extended structure block, extract the extended structure contour and call the slip line segments within the contour. Arrange the line segments in sequence as time progresses and record the displacement direction. According to the displacement direction, combine the blocks corresponding to the line segments into a continuous slip trajectory in time sequence. Call the trajectory group to construct a spatiotemporal overlay graphic and generate the spatiotemporal evolution result of the rock mass by natural collapse method.

[0054] The joint closure configuration block includes the spatial distribution of intersection points, the network of closed path lines, and the configuration boundary outline. The deformation trend correlation block includes the deformation rate difference zone, the direction change trend zone, and the continuously evolving boundary circle. The slip path control block includes the consistent slip node group, the slip direction dominant line, and the node connection sequence line segment. The stable boundary extension structure block includes the path intersection central area, the propagation channel integration area, and the extended configuration boundary. The spatiotemporal evolution results of rock mass by natural collapse method include the time series slip trajectory, the directional displacement annotation group, and the spatiotemporal composite evolution graphics.

[0055] Please see Figure 2 The specific steps of S1 are as follows:

[0056] S101: Based on the spatial orientation of joint surfaces in the three-dimensional structure of rock mass, extract the intersection coordinates of intersecting joint surfaces, calculate the spatial distance and orientation angle between intersection points, and filter them by combining preset angle and distance thresholds to generate a sequence of spatial positions of intersection points;

[0057] First, it is necessary to identify and extract parameters from each joint surface in the 3D rock mass structure model. The spatial attribute information of each joint surface in the model is retrieved one by one, including its planar azimuth, dip angle, center point coordinates, and boundary vertex coordinates. The spatial position data of each joint is recorded face by face through the model data interface. Based on this, the axis-aligned bounding box method is used to determine whether any two joint surfaces have a potential intersection. For joint surfaces with the possibility of intersection, the point set of their boundary contours is densified. A dense point sequence is generated at 0.5-meter intervals for each boundary line for subsequent geometric calculations. Then, a spatial equation system is constructed for adjacent joint surfaces to obtain their theoretical intersection line. Next, boundary conditions are used to determine whether the intersection line actually intersects the two surfaces. If the geometric conditions are met, the two endpoints of the intersection line are retained as the meeting point, and its 3D coordinates in the model space are retrieved and saved again. For each pair of meeting points... The spatial distance between them is calculated and the corresponding horizontal projection angle is recorded. A threshold judgment is made on the distance between each pair of intersection points. The threshold is set to 3 meters. That is, when the distance between two points does not exceed 3 meters, they are considered to have spatial proximity. At the same time, the difference between their horizontal angle and the main direction of the target joint surface is calculated. If the difference does not exceed 15 degrees, the pair of intersection points is retained as a valid point pair. In actual operation, if the intersection point p1 is (10.0, 5.0, 3.0) and the intersection point p2 is (12.5, 6.8, 3.2), the spatial distance is 3.1 meters. If the condition is not met, the point pair is removed. If the intersection point p3 is (10.4, 5.2, 3.1), the distance between it and p1 is 0.6 meters and the direction difference is 12 degrees. If the condition is met, it is retained. Finally, a series of spatial relationships between valid intersection points are constructed in this way and arranged in order to form a sequence of spatial positions of intersection points.

[0058] S102: Call the spatial location sequence of intersection points, calculate the difference between Euclidean distance and direction, remove inconsistent point pairs based on spatial closure threshold and included angle threshold, construct an adjacency matrix that meets the connection conditions, and generate a node map of joint closable paths.

[0059] First, the 3D coordinates of point pairs are extracted from the sequence. The Euclidean distance between each pair is calculated, and this distance is used as one of the criteria for judgment. Simultaneously, the orientation value of the joint surface to which the point pair belongs is retrieved, and the difference between these orientation values ​​is calculated to obtain the orientation difference. If the difference between the two values ​​does not exceed a set angle threshold of 20 degrees, the point pair is considered to have approximately the same orientation. After performing the above calculations on all point pairs, point pairs that do not meet the dual constraints of distance and angle are eliminated. Then, an adjacency matrix is ​​constructed using the remaining point pairs that meet the conditions. Each element in the matrix represents a connectable relationship between two points. If the distance between point pairs is less than or equal to 3 meters and the orientation difference is less than or equal to 20 degrees, then the pair is considered to be approximately connected. If the distance between points p1 and p2 is 2.6 meters and the direction difference is 17 degrees, then the position is marked as 1 in the matrix; otherwise, it is 0. In actual construction, for example, if there is a distance between points p1 and p2 of 2.6 meters and the direction difference is 17 degrees, then this pair of points is retained and A[1][2]=1 in the matrix; if the distance between points p2 and p3 is 3.5 meters and the direction difference is 18 degrees, then it is excluded because it exceeds the distance threshold and A[2][3]=0 in the matrix. After performing complete calculation and judgment on all intersection pairs, the adjacency matrix is ​​constructed. The number of rows and columns of the matrix is ​​equal to the number of intersections. The final adjacency matrix describes the connectable relationship of each intersection under the constraints of spatial direction and distance, and provides input data for the subsequent generation of path map.

[0060] S103: Based on the node map of closable joint paths, extract the closed paths and sort them according to spatial continuity. Connect the node coordinates in the paths, splice them to form a closed boundary profile, and generate joint closure configuration blocks.

[0061] By traversing all nodes in the graph and using a sequential connection method, all possible sequences of points that could form a closed path are found. Starting from any point, the path with a value of 1 in the adjacency matrix is ​​sequentially connected to the next unvisited point, continuing until a closed path returning to the starting point is formed. Throughout the path, the distance between any adjacent nodes must be less than 3 meters, the angular difference between any adjacent joint surfaces must be less than 20 degrees, and the spatial distance between the first and last points must be less than 5 meters to be considered a valid closed path. If these conditions are not met, the path is discarded. For example, starting from point p1 = (10.0, 5.0, 3.0), the path is sequentially connected to p2 = (10.8, 5.6, 3.1) and p3 = (11.4...). (6.2, 3.0), then return to p1 to form a three-point closed path. It is determined that the distance of each segment is within 2.0 meters, the maximum angle difference is 15 degrees, and the distance between the first and last ends is 1.8 meters. All conditions are met, so this path is retained. All paths that meet the conditions form a path set. Next, the nodes in the path are arranged in the order of connection. The coordinates of these points are connected in sequence to form a continuous closed line. The point set between the line segments is supplemented by coordinate interpolation to form a smooth curve. Finally, the closed curve is used as the boundary to generate a three-dimensional joint configuration block. This block is represented as a spatial graphic spliced ​​by several closed paths. Its boundary line corresponds to the intersection point path, thus completing the splicing generation of the configuration block.

[0062] Please see Figure 3 The specific steps of S2 are as follows:

[0063] S201: Based on the joint closure configuration block, extract the deformation monitoring coordinates of the corresponding area, collect the displacement values ​​of the coordinate points at different times, calculate the deformation rate per unit time, and generate the regional deformation rate matrix.

[0064] First, the spatial boundary coordinates of the configuration blocks are retrieved. For each block region, the coverage area of ​​the monitoring points is determined. Deformation monitoring points falling within this boundary area are selected, and the initial three-dimensional coordinates of each point and the corresponding coordinate changes in each subsequent monitoring cycle are recorded sequentially. The time interval is set in days, with common settings such as 5 days or 7 days. For each monitoring point, starting from the first moment, its subsequent coordinate values ​​are read sequentially. The displacement difference in the X, Y, and Z directions between two adjacent moments is calculated. The differences in the three directions are squared, added, and squared to obtain the total displacement within the corresponding monitoring cycle. Then, division is performed by time interval to obtain the deformation rate per unit time within that period. Record the deformation rate of the point in each monitoring cycle. In the example, if point P1 is (102.3, 86.5, 14.2) in the first period and (102.5, 86.8, 14.3) in the second period, with a time interval of 5 days, and the displacement is calculated to be approximately 0.39 meters, then its deformation rate is 0.078 meters / day. If there are 10 monitoring points in a certain area, then the deformation rate of each point in each period is arranged into a 10-row N-column numerical table, where each column represents a certain period and the row represents the corresponding monitoring point. In the final output, it is saved in matrix form, that is, the regional deformation rate matrix is ​​generated. Each value in the matrix corresponds to the deformation rate of a certain monitoring point per unit time at a certain moment.

[0065] S202: Call the regional deformation rate matrix, calculate the rate difference between adjacent points and obtain the direction vector, record the points where the direction difference remains consistent across multiple time series, filter the set of continuously shifting points, and obtain the index set of continuously changing regions.

[0066] First, each row in the matrix is ​​processed. For adjacent data points, the difference in deformation rate is calculated, extracting the rate change of the monitoring point across multiple time periods. Then, each monitoring point is compared pairwise with its adjacent points to determine if the difference in deformation rate between adjacent points within the same time period is significant. A significant rate difference threshold of 0.05 m / day is set; if the rate difference between two points exceeds this value in a given period, a rate difference is considered to exist. In the example, if point A's rate is 0.03 m / day and point B's is 0.09 m / day, the difference is 0.06, meeting the threshold condition, so the point pair is retained. Further analysis of the displacement direction of this point pair across all time periods is performed. The direction vector from point A to point B is recorded for each time period, and its projection angle on the horizontal plane is calculated. For three or more consecutive periods of direction differences... In comparison, the criteria for determining whether the direction difference remains consistent are as follows: if the direction difference between any two adjacent time periods does not exceed 10 degrees, the point pair is considered to have a consistent direction. For example, if the direction of a point pair is 80 degrees in the first period, 84 degrees in the second period, and 85 degrees in the third period, then the direction is considered to be continuous and stable. If the rate difference meets the set requirements and the direction difference is continuous and stable, then the two points in the point pair are marked as points with consistent directions, and these points are saved into a candidate set. Further judgment is made on all points in the candidate set, that is, the deformation rate of the same monitoring point in three consecutive time periods is analyzed for trends. If the rate increases or decreases continuously, then the point is marked as a trend point. For example, if the rate of point C is 0.01, 0.02, and 0.03 meters / day, it is considered a growth trend point. Finally, these points and their corresponding index numbers are output as a set of indexes for continuously changing areas.

[0067] S203: Based on the continuously changing region index set, establish spatial adjacency relationships and construct closed boundaries, connect coordinates to form contour boundaries, and generate deformation trend related blocks;

[0068] First, extract the 3D coordinates of all marked continuously changing points to construct a point distribution map. Then, sequentially determine the spatial relationship between each pair of points, defining the adjacency criteria as a horizontal distance not exceeding 5 meters and a vertical difference not exceeding 2 meters between points. Calculate each pair of points that meets the criteria, mark them as adjacent points, and record the connection relationship. Based on the adjacent point pairs, construct a spatial connection graph, setting all points as nodes and adjacency relationships as edges. Connect all points that meet the criteria to form several paths. Traverse the paths to find if a closed structure exists; that is, if starting from any point, one can return to the starting point after several adjacent jumps, then a closed structure is considered to be formed. Record the coordinates of all nodes in the closed structure according to their connection order. The boundary segments are formed by sequentially connecting them. Then, the coordinates of the boundary segments in the closed path are connected sequentially to form a complete contour boundary. The boundary segments are smoothed by using sparse point completion to make the boundary contour a continuous graphic. For example, starting from point P1=(100.1, 50.2, 10.5), it is connected sequentially to P2=(101.3, 50.8, 10.6), P3=(101.5, 51.9, 10.5), P4=(100.2, 51.7, 10.3), and finally back to P1 to form a closed region. After generating the contour boundary, the boundary is recorded in the block data structure and assigned a block number. A correspondence is established with the corresponding joint configuration block. Finally, the block is output as a deformation trend association block.

[0069] Please see Figure 4 The specific steps of S3 are as follows:

[0070] S301: Based on the deformation trend association map, obtain the node position of the shear surface in the area, compare the angle between the node sliding direction and the contact direction, and filter the node group with the same direction according to the direction difference threshold. Arrange the nodes into a linear sequence according to the connection order to generate a node sequence with the same direction.

[0071] First, select the shear plane information corresponding to the area covered by the tile. Extract the coordinate information of all nodes related to shearing from the spatial data, obtain the 3D position of each node and the corresponding sliding direction vector. Then, calculate the angle difference between the contact surface direction vector and the sliding direction. For each node, calculate the angle between its sliding direction and the normal vector of the shear plane. Set the direction difference threshold to 15 degrees. When the angle is less than or equal to this threshold, the node is considered to have the same sliding direction as the contact direction, and such nodes are added to the candidate set. In practice, if node A has a sliding direction of 82 degrees and a contact direction of 87 degrees, the direction difference is... If the angle is 5 degrees, less than the threshold, the node is retained. If node B has a sliding direction of 75 degrees, a contact direction of 97 degrees, and an included angle of 22 degrees, which is greater than the threshold, it is removed. After the filtering is completed, the spatial positions of the retained nodes are sorted according to the direction of the shear plane. The nodes are sorted by reading the difference between the horizontal and vertical coordinates, and are first sorted in ascending order in the X-axis direction. If the X values ​​are the same, they are sorted in ascending order in the Y-value to form a continuous linear structure. After sorting all the nodes that meet the filtering requirements, the nodes are connected one by one in a one-dimensional linear sequence according to their spatial connection relationship. The node number, position, sliding direction and direction difference are recorded. The final output is a sequence of nodes with consistent directions.

[0072] S302: Call the sequence of nodes with the same direction, calculate the direction vector of adjacent nodes, and determine whether they constitute a valid line segment based on the sliding advance threshold. The line segments that meet the advance conditions are grouped into a set according to the node order to obtain the sliding advance line segment set.

[0073] Starting from the first node, the coordinates between it and the next node are read in pairs. The direction vector between each pair of adjacent nodes is calculated, and the continuity of the sliding direction is determined by the difference. The start point, end point, and direction angle of each vector segment are recorded. A threshold of 12 degrees is set for the angle difference in sliding propulsion; that is, if the angle difference between two consecutive direction vectors is less than or equal to 12 degrees, they are considered to constitute a sliding propulsion relationship. Furthermore, it is necessary to determine whether the direction of the line segment is close to the original main sliding direction. The main direction angle is set as the target θ0. If the difference between the direction angle and θ0 does not exceed 15 degrees, the propulsion direction is considered consistent. In the example... If the direction from node P1 to P2 is 78 degrees, and the direction from P2 to P3 is 82 degrees, with a direction difference of 4 degrees, and the target direction θ0 is 80 degrees, with a deviation of 2 degrees from this direction, then this line segment is considered valid, and the next segment is added. If the direction difference of a certain segment exceeds the set threshold, such as P3 to P4 being 96 degrees, with a direction difference of 14 degrees exceeding 12 degrees, then the advancement is interrupted, the current line segment combination ends, all line segments that meet the conditions are summarized, and their starting and ending coordinates, direction angles, included angle differences, and corresponding node index numbers are recorded. All line segments that meet the advancement conditions are combined into a line segment set according to their original node connection order, and finally, the sliding advancement line segment set is obtained.

[0074] S303: Based on the set of sliding advance line segments, retrieve the line segment direction vectors and aggregate them according to the consistency of direction, gather continuous line segments into a single path structure, and generate the sliding path master control block;

[0075] The direction vectors of all line segments are extracted one by one, and a direction consistency judgment is performed. First, the line segments are sorted according to their starting coordinates. Adjacent line segments are combined one by one according to the principle of spatial distance priority. It is judged whether the angle difference between the direction of the current line segment and the direction of the next line segment is lower than the set direction consistency threshold, which is set to 10 degrees. If the condition is met, the two line segments are merged into the same structural path, and the starting point, ending point and direction attributes of the merged path are recorded. In the continuous judgment process, if the direction difference between three consecutive line segments all meet the consistency condition, it is identified as a continuous main control path and added to the merging sequence. In actual processing, if the line segments If L1 is 80 degrees, L2 is 83 degrees, and L3 is 87 degrees, with directional differences of 3 and 4 degrees respectively, all less than the threshold, then the three segments can be continuously aggregated. If L4 is 96 degrees, with a difference of 9 degrees from L3, which is also less than the threshold, then aggregation continues. If L5 is 105 degrees, with a difference of 9 degrees from L4, but the cumulative deviation from L1 exceeds 25 degrees, then aggregation stops, forming a main path structure block. For all continuously aggregated line segment combinations, record their total length, average direction angle, number of nodes, and path number. Finally, express these continuous structure paths in three-dimensional space and output them in the form of blocks, that is, generate the sliding path master control block.

[0076] Please see Figure 5 The specific steps of S4 are as follows:

[0077] S401: Based on the sliding path master control block, extract the intersection coordinates between all path lines in the path line set, count the distribution frequency of the intersection coordinates in the block, and construct a spatial distribution heat map based on the density of coordinates. Extract the contours of continuous clustered areas in the heat map to generate the spatial distribution range of the intersection area.

[0078] First, the coordinate data of all path lines in the tile is retrieved, and possible intersection points between the path lines are extracted. It is then determined whether a spatial intersection exists between each pair of path lines; if so, its intersection coordinates are extracted, and all intersection coordinates are uniformly incorporated into a 3D spatial indexing system. The spatial position of each intersection point is recorded in the form of (x, y, z). The frequency of occurrence of identical or adjacent coordinate points is obtained by counting the number of times they appear. During the statistical process, the entire tile space is divided into equally sized cubic grid units, each with a side length of 1 meter. The number of intersection points in each grid is counted as a heat index. When a grid contains an intersection point... When the number is ≥3, the grid is considered a hotspot region. The hotspot grids are further aggregated into continuous regions according to their positions. The contour extraction operation is performed on the region formed by connecting multiple consecutive hotspot units. The edge points are extracted using the grid boundary connection principle. By identifying the connection relationship between each hotspot unit and its adjacent units, a set of contour points in the heat map is generated. In the example, if the total number of intersection points in a certain block is 420, there are 18 hotspot grids after division. The number of intersection points of each hotspot grid is between 3 and 9. These hotspot grids are then connected in space to generate a continuous heat distribution region, and a closed contour is formed according to its outer boundary coordinates. The final output is the spatial distribution range of the intersection area.

[0079] S402: Call the spatial distribution range of the intersection area, determine the number of overlapping paths and the density of intersection points within the intersection area, mark the area with an overlap rate exceeding the path overlap threshold as a concentrated propagation area, and extract the boundary by combining the path extension direction to obtain the boundary of the concentrated propagation area of ​​the path.

[0080] The spatial trajectories of all path lines within the specified range are traversed, determining whether each path line crosses an intersection area. Within the intersection area, each path line is segmented and its segments are counted, calculating the number of segments that the path line crosses. The number of overlapping paths in this area is also recorded. If the number of overlapping paths in an area exceeds a set path overlap threshold, it is identified as a concentrated propagation area. Here, the overlap threshold is set to 30% of the total number of paths. If the total number of path lines in the map is 40, the threshold is 12. If 14 path lines cross the center of an area, the area is considered to meet the concentrated propagation judgment condition. After the judgment is completed... Record the numbers of all intersection areas that meet the overlap requirements in sequence. Then, extract the boundary of the region based on the direction information of the path lines. The extension direction is determined based on the value of the main direction vector of each path line. The direction difference tolerance is set to 15 degrees. That is, if the direction consistency of multiple paths at a certain boundary meets the condition that the difference does not exceed 15 degrees, then the direction is taken as the extension direction. The boundary is extended along the direction. The extension length is set to 50% of the average diameter of the intersection area. For example, if the average diameter of an intersection area is 20 meters, then the extension length is 10 meters. During this extension process, record the coordinates of the newly formed boundary points to complete the generation of the boundary of the path concentration propagation area.

[0081] S403: Based on the boundary of the path-concentrated propagation area, aggregate the coordinates of control points on the outline of the area, arrange the directions and reconstruct the node connections according to the principle of consistency of direction, form a closed structural framework and perform boundary fitting to generate stable boundary extension structural blocks;

[0082] The coordinates of control points on the boundary line are extracted point by point and summarized into a coordinate point set. For each control point, its sequential number and spatial position parameters on the boundary curve are recorded. The spatial vector between each pair of adjacent control points is retrieved, and their direction angle values ​​are calculated. All direction angles are sorted sequentially to determine consistency. If the difference between three or more consecutive directions is within 10 degrees, the sequence is considered a consistent direction segment. Nodes are then connected and reconstructed in this order. During the connection and reconstruction process, points that are too close are merged. If the distance between two points is less than 0.3 meters, they are merged into one control point; if the distance is between 0.3 meters and 1.0 meter... If a node is retained as a single node, a continuous chain of nodes with consistent direction is formed after processing all control points. The chain is then closed. If the distance between the first and last nodes is less than 3 meters, the chain is automatically closed. If it exceeds 3 meters, interpolation is performed to generate a closed loop structure. The structure is then fitted to its boundary. Piecewise linear fitting is used to record the start and end points and direction angles of each segment. The fitting accuracy is controlled within ±1 meter error range. For example, if a path contains 26 control points in its concentrated propagation boundary, 18 of them are retained as reconstruction reference points after screening. After fitting, a closed structural framework is constructed, and the framework is output as a stable boundary extension structure block.

[0083] Please see Figure 6 The specific steps of S5 are as follows:

[0084] S501: Based on the stable boundary extended structure block, extract the extended structure outline within the block, call the time series labels of the sliding line segments within the outline, rearrange the line segments in the time progression order, and record the spatial displacement direction corresponding to the line segments to generate a time series sliding line segment set;

[0085] First, extract all slip structure information contained within the closed boundaries of the tile. Identify and read the set of contour lines constituting the extended structure, and obtain the slip segments contained within each contour line. Sequentially retrieve the unique identifier, endpoint coordinates, and associated time stamp information of these segments. This stamp is usually a date stamp or a unique number. Sort the segments in ascending order of timestamps from earliest to latest. Reorder all segments according to their time stamp values. After sorting, extract the slip direction data for each segment in sequence. Record the vector difference between endpoints as an angle. Each segment must be labeled. The corresponding sliding direction angle is calculated. For example, if line segment L1 slides from P1 (102.3, 85.7) to P2 (104.8, 87.9), the direction angle is calculated to be 56 degrees. Then, the direction of L1 is recorded as 56 degrees and it is placed first in the sequence according to time. If the time label of L2 is later than that of L1 and the direction is 58 degrees, it is placed after L1 and the sequence continues to be built. After all the sliding line segments are arranged according to time, they are organized into a linear structure. The start and end coordinates, time label and direction angle of each line segment are recorded. Finally, a time series sliding line segment set composed of time progression order is generated.

[0086] S502: Call the time series sliding line segment set, merge the corresponding blocks of the line segments within the same contour in time order, and maintain the consistency constraint on the displacement direction of each line segment. Connect the sequence blocks to construct a continuous sliding path and obtain the sliding trajectory sequence structure.

[0087] All line segments within the same contour boundary are combined according to their time sequence. Starting with the earliest line segment, subsequent time nodes' line segments are sequentially pieced together to construct a path block. For each line segment, a direction consistency check is performed. The direction angle value between the current line segment and its preceding line segment is extracted, the angle difference is calculated, and compared with a set direction consistency threshold. If the angle difference is less than or equal to 10 degrees, the direction is considered consistent, and the segments can be connected into a continuous path. If the difference exceeds 10 degrees, the connection is interrupted, and the segments are recorded separately. For example, if line segment L1 has a direction of 57 degrees and L2 has a direction of 63 degrees, the difference is 6 degrees, satisfying the condition, and the segment is connected. If L3 has a direction of 74 degrees, the difference with L2 is... If the angle exceeds 11 degrees, a new path segment is started. This logic is used to traverse the entire sequence set. For each set of validly spliced ​​line segments, a continuous structure block is generated. The start and end positions, sliding direction, and time labels of each line segment are recorded synchronously during the splicing process and arranged in chronological order to form a complete path sequence. The coordinates of the connecting parts between the blocks are aligned. If there is a connection offset greater than 0.5 meters, a midpoint calibration operation is performed. The starting position of the subsequent line segments is adjusted according to the starting point to control the splicing error within a reasonable range. Finally, multiple continuous sliding paths are formed. Each path consists of multiple sliding line segments, forming a structured sliding trajectory sequence structure.

[0088] S503: Based on the slip trajectory sequence structure, extract the spatial boundaries and time labels of trajectory blocks, establish the spatial superposition relationship of trajectories with different time periods, and generate the spatiotemporal evolution results of rock mass using the natural collapse method.

[0089] First, read the boundary coordinate data and time label information of all trajectory path tiles. Group the path sequences by time label, with each group corresponding to the trajectory status within a time period. Perform overlay processing on the path tiles of each time period according to the spatial boundaries, extracting the outermost coordinate set of the path boundaries within each time period. Compare the boundary changes between the current time period and the previous time period according to the overlay order to determine whether the trajectory has undergone spatial expansion, contraction, or turning. For spatially expanded areas, record the newly added boundary coordinates incrementally. For overlapping areas, determine whether the number of paths has changed. If the number increases and the direction is the same, mark it as repeated advancement; if the direction is different... Conversely, it is marked as a slip conflict area. In the example, if the trajectory boundary of time period T1 is (100, 90)-(120, 110) and T2 is (102, 92)-(124, 115), it is judged that the boundary is expanding in the southeast direction. The expansion trend and the amount of change are recorded. Data of multiple time periods are processed in this way. The change status between the boundaries of adjacent time periods is compared to establish the spatial superposition relationship of the trajectory of each stage in the time series. All superposition results are encoded and output. Time labels and spatial coverage area identifiers are added to each path segment. Finally, the spatiotemporal evolution results of the rock mass of the natural collapse method are output based on this superposition structure.

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

Claims

1. A method for simulating the spatiotemporal evolution of rock masses using the natural collapse method, characterized in that, Includes the following steps: S1: Extract the orientation of the three-dimensional joint surfaces of the rock mass, obtain the location of the intersection points, compare the spatial distance and direction between the intersection points, filter the closable path lines, sort them according to spatial continuity and splice the paths to generate joint closure configuration blocks; S2: Based on the joint closure configuration block, extract the regional deformation rate, compare the rate difference and direction of adjacent regions, filter regions with continuous difference direction, and generate deformation trend association blocks; S3: Based on the deformation trend association block, extract the shear surface node, compare the difference between the node slip and contact direction, filter the node group with the same direction and sort them into line segments, determine the slip propulsion direction, and generate the slip path main control block; S4: Extract the path intersection positions based on the sliding path master control block, count the number of intersections and verify the intersection area range, determine the concentrated propagation area, reorganize to form an extended structure outline, and generate a stable boundary extended structure block; S5: Extract slip segments from the stable boundary extension structure block, arrange them in time progression order and record the displacement direction, and combine the line segment blocks into slip trajectory according to time sequence to generate the spatiotemporal evolution result of rock mass by natural collapse method; The region where the difference direction is continuous refers to the region where the direction of the difference in deformation rate between adjacent regions is consistent and continuously distributed in space. The specific steps of S1 are as follows: S101: Based on the spatial orientation of joint surfaces in the three-dimensional structure of rock mass, extract the intersection coordinates of intersecting joint surfaces, calculate the spatial distance and orientation angle between intersection points, and filter them by combining preset angle and distance thresholds to generate a sequence of spatial positions of intersection points; S102: Call the spatial location sequence of the intersection point, calculate the difference between the Euclidean distance and the direction, remove inconsistent point pairs according to the spatial closure threshold and the included angle threshold, construct an adjacency matrix that meets the connection conditions, and generate a node map of joint closable paths. S103: Based on the node map of the joint closable path, extract the closed path and sort it according to spatial continuity, connect the node coordinates in the path, splice them to form a closed boundary contour, and generate a joint closure configuration block.

2. The method for simulating the spatio-temporal evolution of a rock mass subjected to natural caving according to claim 1, characterized in that, The joint closure configuration block includes the spatial distribution of intersection points, the network of closed path lines, and the configuration boundary outline. The deformation trend association block includes the deformation rate difference zone, the direction change trend zone, and the continuously evolving boundary circle. The slip path master control block includes the consistent slip node group, the slip direction dominant line, and the node connection sequence line segment. The stable boundary extension structure block includes the path intersection central area, the propagation channel integration area, and the extended configuration boundary. The spatiotemporal evolution results of the rock mass obtained by the natural collapse method include the time series slip trajectory, the directional displacement annotation group, and the spatiotemporal composite evolution graph.

3. The method according to claim 1, characterized in that, The concentrated propagation area refers to the area where slip paths frequently intersect in space and slip activities occur in a highly concentrated manner.

4. The method for simulating the spatiotemporal evolution of rock masses using the natural collapse method according to claim 1, characterized in that, The specific steps of S2 are as follows: S201: Based on the joint closure configuration block, extract the deformation monitoring coordinates of the corresponding region, collect the displacement values ​​of the coordinate points at different times, calculate the deformation rate per unit time, and generate a regional deformation rate matrix. S202: Call the region deformation rate matrix, calculate the rate difference between adjacent points and obtain the direction vector, record the points where the direction difference remains consistent in multiple time series, filter the set of continuously shifting points, and obtain the index set of continuously changing regions. S203: Based on the continuously changing region index set, establish spatial adjacency relationships and construct closed boundaries, connect coordinates to form contour boundaries, and generate deformation trend association blocks.

5. The method for simulating the spatio-temporal evolution of a rock mass subjected to natural caving according to claim 1, characterized in that, The specific steps for S3 are as follows: S301: Based on the deformation trend association block, obtain the node position of the shear surface in the area, compare the angle between the node sliding direction and the contact direction, and filter the node group with the same direction according to the direction difference threshold. Arrange the nodes into a linear sequence according to the connection order to generate a node sequence with the same direction. S302: Call the sequence of nodes with the same direction, calculate the direction vector of adjacent nodes, and determine whether they constitute a valid line segment based on the sliding advance threshold. The line segments that meet the advance conditions are grouped into a set according to the node order to obtain the sliding advance line segment set. S303: Based on the set of sliding propulsion line segments, retrieve the line segment direction vectors and aggregate them according to the consistency of direction, gather the continuous line segments into a single path structure, and generate the sliding path master control block.

6. The method for simulating the spatio-temporal evolution of a rock mass subjected to natural caving according to claim 1, characterized in that, The specific steps of S4 are as follows: S401: Based on the sliding path master control block, extract the intersection coordinates between all path lines in the path line set, count the distribution frequency of the intersection coordinates in the block, and construct a spatial distribution heat map based on the density of coordinates. Extract the contours of continuous clustered areas in the heat map to generate the spatial distribution range of the intersection area. S402: Call the spatial distribution range of the intersection area, determine the number of overlapping paths and the density of intersection points within the intersection area, mark the area with an overlap rate exceeding the path overlap threshold as a concentrated propagation area, and extract the boundary in combination with the path extension direction to obtain the boundary of the concentrated propagation area of ​​the path. S403: Based on the boundary of the path-centralized propagation area, the coordinates of control points on the contour line of the aggregated area are used to arrange the directions and reconstruct the node connections according to the principle of consistency of direction, forming a closed structural framework and performing boundary fitting to generate stable boundary extension structural blocks.

7. The method for simulating the spatiotemporal evolution of rock masses using the natural collapse method according to claim 1, characterized in that, The specific steps of S5 are as follows: S501: Based on the stable boundary extended structure block, extract the extended structure outline within the block, call the time series labels of the sliding line segments within the outline, rearrange the line segments according to the time progression order, and record the spatial displacement direction corresponding to the line segments to generate a time series sliding line segment set; S502: Call the time series sliding line segment set, merge the corresponding blocks of the line segments in the same contour in time order, and maintain the consistency constraint on the displacement direction of each line segment. Connect the sequence blocks to construct a continuous sliding path to obtain the sliding trajectory sequence structure. S503: Based on the slip trajectory sequence structure, extract the spatial boundaries and time labels of the trajectory blocks, establish the spatial superposition relationship of the trajectories in different time periods, and generate the spatiotemporal evolution results of the rock mass by the natural collapse method.

8. The method according to claim 7, characterized in that, The continuous sliding path refers to a sequence of sliding line segments that are continuous in time, have the same displacement direction in space, and are formed by piecing together blocks.

Images