Potential safety hazard analysis method and system based on underground cavern three-dimensional geological intelligent modeling

Abstract

The invention provides a potential safety hazard analysis method and system based on underground cavern three-dimensional geological intelligent modeling, and the method comprises the steps: firstly, based on underground cavern three-dimensional model data, cutting the underground cavern three-dimensional model data into a plurality of sub-models according to a geological condition segmentation rule, and enabling the geological condition segmentation rule to dynamically adjust the cutting spacing according to the gradient change characteristics of a rock stratum; performing gradient analysis on the longitudinal and axial rock stratum attributes of each sub-model to generate rock stratum transition characteristic data containing the adjacent rock stratum attribute change rate and the transition area space distribution, and dynamically identifying the surrounding rock type of the sub-model according to the rock stratum transition characteristic data; and mapping and matching the surrounding rock identification result with a preset protected object space coordinate, collecting a sub-model geological parameter set in real time, finally triggering a safety early warning signal according to a geological parameter set numerical value interval, and generating an emergency disposal scheme set associated with the safety early warning signal, so that the potential safety hazard of the underground cavern can be effectively analyzed.

Description

Technical Field

[0001] The present invention relates to the field of computer technology, and more specifically, to a method and system for analyzing potential safety hazards based on three-dimensional geological intelligent modeling of underground caverns. Background Art

[0002] In the field of underground cavern engineering, with the continuous expansion of project scale and the increasing complexity of geological conditions, the accurate analysis and effective response to potential safety hazards have become a crucial link in ensuring project safety. However, most of the existing methods for analyzing potential safety hazards rely on traditional geological exploration means and empirical judgments, and have many limitations.

[0003] On the one hand, the existing three-dimensional geological modeling technologies often adopt a fixed grid division method and cannot be dynamically adjusted according to the actual changes in geological conditions. This static modeling method is difficult to accurately reflect the complex changes in the rock formations around underground caverns, especially in areas with large rock formation gradient changes, which easily leads to insufficient model accuracy and thus affects the accuracy of potential safety hazard analysis.

[0004] On the other hand, the existing methods for potential safety hazard warning and disposal often lack real-time association with the set of geological parameters and are difficult to dynamically adjust the warning threshold and disposal measures according to the real-time changes in geological conditions. This static warning and disposal mechanism often proves inadequate in the face of complex and changeable geological environments and cannot effectively respond to potential safety hazards in a timely manner. Summary of the Invention

[0005] In view of the above-mentioned problems, in combination with the first aspect of the present invention, embodiments of the present invention provide a method for analyzing potential safety hazards based on three-dimensional geological intelligent modeling of underground caverns, the method including: Based on the three-dimensional model data of the underground cavern, cutting the three-dimensional model into multiple sub-models according to the geological condition segmentation rule, and the segmentation rule dynamically adjusts the cutting spacing of each sub-model according to the rock formation gradient change characteristics; Performing gradient analysis on the longitudinal rock formation attributes and axial rock formation attributes of each sub-model to generate rock formation transition characteristic data of the sub-model, where the rock formation transition characteristic data includes the attribute change rate between adjacent rock formations and the spatial distribution of the transition area; According to the rock formation transition characteristic data, dynamically identifying the surrounding rock type of each sub-model, and mapping and matching the surrounding rock identification result with the spatial coordinates of a preset protection object to obtain a mapping and matching relationship, and based on the mapping and matching relationship, collecting the set of geological parameters of the sub-model in real time; Triggering a safety warning signal according to the numerical range of the set of geological parameters, and generating a set of emergency disposal plans associated with the safety warning signal.

[0006] In another aspect, an embodiment of the present invention further provides a safety hazard analysis system based on three-dimensional geological intelligent modeling of underground chambers, including a processor and a machine-readable storage medium. The machine-readable storage medium is connected to the processor and is used to store programs, instructions or codes. The processor is used to execute the programs, instructions or codes in the machine-readable storage medium to implement the above method.

[0007] Based on the above aspects, the embodiment of the present invention realizes the efficient and accurate analysis of safety hazards in underground chambers, significantly improving the safety and reliability of underground chamber projects. Specifically, first, according to the three-dimensional model data of the underground chamber and in combination with the dynamically adjusted geological condition segmentation rules, the complex three-dimensional geological model is finely cut into multiple sub-models, which not only considers the gradient change characteristics of the rock strata but also ensures that each sub-model can accurately reflect the subtle differences in local geological conditions. On this basis, further gradient analysis is carried out on the longitudinal rock stratum attributes and axial rock stratum attributes of each sub-model to generate rock stratum transition characteristic data including the attribute change rate between adjacent rock strata and the spatial distribution of the transition area, which not only reveals the internal relationship and change law between the rock strata but also provides a key basis for the dynamic identification of surrounding rock types. By mapping and matching the surrounding rock identification results with the spatial coordinates of the preset protected object, a set of geological parameters closely related to the protected object can be collected in real time, thereby realizing the accurate positioning and real-time monitoring of safety hazards. Finally, according to the numerical range of the set of geological parameters, the method can intelligently trigger a safety warning signal and generate an associated set of emergency response plans, improving the response speed to safety hazards. In summary, a comprehensive and in-depth analysis and efficient response to safety hazards in underground chambers are realized, which helps the safe construction and operation of underground chamber projects. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 is a schematic flowchart of the execution process of the safety hazard analysis method based on three-dimensional geological intelligent modeling of underground chambers provided by an embodiment of the present invention.

[0009] Figure 2 is a schematic diagram of exemplary hardware and software components of the safety hazard analysis system based on three-dimensional geological intelligent modeling of underground chambers provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0010] The present invention will be specifically described below in conjunction with the accompanying drawings of the specification. Figure 1 is a flowchart of the safety hazard analysis method based on three-dimensional geological intelligent modeling of underground chambers provided by an embodiment of the present invention. The safety hazard analysis method based on three-dimensional geological intelligent modeling of underground chambers will be introduced in detail below.

[0011] Step S110: Based on the 3D model data of the underground cavern, cut the 3D model into multiple sub-models according to the geological condition segmentation rules, where the segmentation rules dynamically adjust the cutting spacing of each sub-model according to the characteristics of the rock layer gradient change.

[0012] For example, in an underground construction project, in order to ensure construction safety and the stability of subsequent operation, it is necessary to model the underground caverns involved and analyze potential safety hazards. Since the geological conditions at the location of the underground construction project are relatively complex, there are multiple different types of rock layers, and the distribution and characteristics of the rock layers vary greatly. Therefore, it is first necessary to obtain the 3D model data of the underground caverns and divide them reasonably.

[0013] Step S111: Generate an initial 3D geological model by fusing geological radar detection data, core drilling data, and UAV oblique photography data. Among them, the geological radar detection data analyzes the distribution of rock mass fractures through the reflection signals of multi-frequency electromagnetic waves, the core drilling data calibrates the rock layer attributes through the physical and mechanical parameters of the core samples, and the UAV oblique photography data reconstructs the surface and shallow geological structures through multi-view images.

[0014] In this embodiment, a geological radar can be used to conduct a comprehensive detection of the planned area of the underground cavern. The geological radar can move on the ground at a set spacing and route, emitting electromagnetic waves of different frequencies. Among them, low-frequency electromagnetic waves (such as 100 MHz) have strong penetration ability and can penetrate deep into the ground, and can be used to detect the general structure of deep rock masses and possible large fractures; high-frequency electromagnetic waves (such as 1000 MHz) have higher resolution but shallower penetration depth, and are mainly used to detect tiny fractures in shallow rock masses.

[0015] When the electromagnetic wave encounters a fracture in the rock mass, reflection will occur, and the characteristics such as the intensity, phase, and propagation time of the reflection signal will change. For example, if the intensity of the reflection signal suddenly weakens and the propagation time prolongs, it may indicate encountering a large fracture; if the reflection signal shows high-frequency oscillation and large intensity fluctuations, it may mean the existence of tiny fractures. By analyzing and processing a large number of reflection signals, the distribution image of the rock mass fractures can be drawn, and information such as the location, trend, and approximate size of the fractures can be determined.

[0016] At the same time, core drilling work can be carried out at different positions in the planned area of the underground cavern. For example, according to the preliminary detection results of the geological radar and the geological characteristics of the area, 20 key positions were selected for drilling. The drilling depth depends on the expected depth of the cavern and the geological conditions, and the deepest reached 50 meters. After taking out the core samples, a series of physical and mechanical parameter tests are carried out on them.

[0017] For the compressive strength test, the core samples are processed into cylinders of standard size and placed in a compression testing machine. Pressure is applied at a constant loading rate until the core samples fail. The maximum pressure value at failure is recorded, and then the compressive strength is calculated based on the cross-sectional area of the core samples. For example, if the cross-sectional area of a certain core sample is 100 square centimeters and the maximum pressure at failure is 500,000 Newtons, then its compressive strength is 500,000 Newtons divided by 100 square centimeters (which needs to be converted to square meters), that is, 50 MPa.

[0018] For the elastic modulus test, the core samples are loaded in stages on a compression testing machine, and the deformation of the samples at different loading stages is measured. The elastic modulus is calculated through the stress-strain curve. The Poisson's ratio test measures the lateral strain while measuring the axial strain, and the ratio of the two is the Poisson's ratio. Through the above tests, the rock formation properties at each drilling position can be accurately calibrated.

[0019] In addition, an unmanned aerial vehicle (UAV) is used for oblique photography. The UAV can fly over the planned area of the underground chamber according to a preset flight path. For example, the flight altitude is set to 50 meters and the flight speed is 5 meters per second. The UAV is equipped with multiple cameras at different angles and takes pictures of the ground from multiple perspectives such as vertical and oblique, obtaining a total of multiple high-resolution images. Through the image matching algorithm, the same feature points in different images can be found and the corresponding relationships between them can be established. Then, using the 3D reconstruction algorithm, based on the corresponding relationships of the above feature points and the internal and external parameters of the cameras, a 3D model of the surface and shallow geological structure is reconstructed. Finally, the data of the rock mass fracture distribution detected by the ground penetrating radar, the rock formation property data obtained from borehole coring, and the data of the surface and shallow geological structure reconstructed by the UAV oblique photography are fused to generate an initial 3D geological model. During the fusion process, according to the spatial coordinate information of the data, the data from different sources can be accurately mapped to the corresponding positions in the 3D space.

[0020] Step S112: Perform spatial grid division on the initial 3D geological model, and based on the rock mass strength distribution curve and the change rate of the permeability coefficient of the grid cells, screen out the non-uniform areas where the difference in rock mass strength and the difference in permeability coefficient between adjacent grid cells exceed the preset thresholds.

[0021] In this embodiment, after generating the initial 3D geological model, it is necessary to perform spatial grid division on it. Specifically, according to the scale and geological characteristics of the underground chamber, the entire initial 3D geological model is divided into grid cells with a specification of 3 meters in length, width, and height. The advantage of this division is that it can not only ensure a detailed description of the geological situation but also prevent the data volume from being too large to cause difficulties in processing.

[0022] For each grid cell, the rock mass strength and permeability coefficient of the grid cell can be calculated using an interpolation algorithm based on the rock formation attribute data obtained from previous drill core sampling and the relevant information detected by ground penetrating radar. The interpolation algorithm takes into account the actual data at the drilling positions around the grid cell and their distance relationships with the grid cell to more accurately estimate the attribute values of the grid cell.

[0023] For example, for a grid cell located in the middle of multiple drill holes, the estimated value of the rock mass strength of the grid cell can be obtained by performing a weighted average according to the inverse of the distance based on the rock mass strength data of the surrounding drill holes. A similar method is used to estimate the permeability coefficient.

[0024] Then, draw the rock mass strength distribution curve for each grid cell. Using the spatial position of the grid cell as the abscissa and the rock mass strength as the ordinate, connect the rock mass strength values of adjacent grid cells to form the rock mass strength distribution curve. At the same time, calculate the change rate of the permeability coefficient between adjacent grid cells, that is, the difference in the permeability coefficients of adjacent grid cells divided by the smaller permeability coefficient value.

[0025] Suppose the preset threshold for the difference in rock mass strength is 12 MPa, and the threshold for the difference in permeability coefficient is 0.015 m / d. Compare all adjacent grid cells one by one. If the difference in rock mass strength between two adjacent grid cells is greater than 12 MPa, or the change rate of the permeability coefficient exceeds 0.015 m / d, mark the area where these two grid cells are located as a non-uniform area.

[0026] For example, the rock mass strength of grid cell A is 75 MPa, and the rock mass strength of its adjacent grid cell B is 90 MPa. The difference between the two is 15 MPa, exceeding the preset threshold of 12 MPa. Then, the area where grid cells A and B are located is marked as a non-uniform area. Similarly, if the permeability coefficient of grid cell C is 0.02 m / d and the permeability coefficient of the adjacent grid cell D is 0.025 m / d, the change rate of the permeability coefficient is (0.025 - 0.02) ÷ 0.02 = 0.25, exceeding the preset threshold of 0.015 m / d. The area where grid cells C and D are located will also be marked as a non-uniform area.

[0027] Step S113: Perform secondary data collection on the non-uniform area to obtain supplementary exploration data on the occurrence of rock joints and the occurrence state of groundwater, and perform point cloud modeling on the excavation surface using a three-dimensional laser scanner. Superimpose the supplementary exploration data and the point cloud modeling data in space with the initial three-dimensional model to generate corrected three-dimensional model data.

[0028] In this embodiment, for the marked non-uniform regions, secondary data collection is required. Specifically, more monitoring devices and measurement points can be added in the above regions to obtain more detailed supplementary exploration data on the occurrence of rock joints and the occurrence state of groundwater.

[0029] When obtaining the data on the occurrence of rock joints, a joint compass is used to measure the rock joints in the non-uniform region in detail. The joint compass can accurately measure the strike, dip, and dip angle of the joints. During the measurement, place the joint compass on the joint surface and read the data according to the correct operation method. For example, at a certain measurement point, the strike of a joint is measured as 30° east of north, the dip is 120° southeast, and the dip angle is 60°. Measure multiple joints in the non-uniform region and record the occurrence data of all joints.

[0030] To obtain the supplementary exploration data on the occurrence state of groundwater, water level observation wells and piezometers can be arranged in the non-uniform region. The water level observation wells are distributed at a certain interval, and by measuring the water level height in the wells, the change of the groundwater level can be understood. The piezometers are installed in the rock mass at different depths to measure the groundwater pressure. For example, in a certain water level observation well, it is measured that the groundwater level is 10 meters from the ground surface; at a certain piezometer, the measured groundwater pressure is 0.1 MPa.

[0031] At the same time, use a three-dimensional laser scanner to scan the excavation surface of the non-uniform region. The three-dimensional laser scanner emits laser beams and quickly obtains the three-dimensional coordinate information of a large number of points on the surface of the excavation surface to form point cloud data. During the scanning, the scanner moves at a certain angle and interval to ensure that the entire excavation surface can be covered.

[0032] Then, process the obtained point cloud data. First, remove the noise points, which may be generated due to external interference or scanner errors. Therefore, a filtering algorithm, such as Gaussian filtering, can be used to smooth the point cloud data and remove the outliers. Then, perform data registration to splice the point cloud data scanned at different positions so that they are accurately aligned in the same coordinate system.

[0033] Finally, spatially superimpose the obtained supplementary exploration data on the occurrence of rock joints, the occurrence state of groundwater, and the point cloud modeling data of the excavation surface with the initial three-dimensional model. During the superimposition process, according to the spatial coordinate information of the data, accurately add the supplementary exploration data and the point cloud modeling data to the corresponding non-uniform regions in the initial three-dimensional model. For example, add the data on the occurrence of rock joints to the model in the form of three-dimensional vectors, display the data on the occurrence state of groundwater in the form of contour lines in the initial three-dimensional model, and fuse the point cloud model of the excavation surface with the model at the corresponding position in the initial three-dimensional model to finally generate the corrected three-dimensional model data.

[0034] Step S114: According to the characteristics of the rock formation gradient change in the corrected 3D model data, cut the sub-model along the axis direction of the cavern with a dynamic spacing. The adjustment logic of the cutting spacing is as follows: when it is detected that the change rate of the rock mass strength in adjacent areas exceeds the set gradient threshold, reduce the cutting spacing to improve the analytical accuracy of the geological attributes of the sub-model, and the included angle between the boundary of each sub-model and the geological structure line is maintained within a set range through spatial vector calculation to avoid acute-angle intersection between the cutting surface and the rock formation strike.

[0035] In this embodiment, after obtaining the corrected 3D model data, the sub-model can be cut along the axis direction of the cavern. First, analyze the characteristics of the rock formation gradient change in the corrected 3D model data, mainly focusing on the change of the rock mass strength, and calculate the change rate of the rock mass strength in adjacent areas, that is, the difference in the rock mass strength between adjacent areas divided by the distance.

[0036] Among them, the gradient threshold can be set to 5 MPa / m. During the cutting process, the change rate of the rock mass strength in adjacent areas is monitored in real time. For example, along the axis direction of the cavern, a monitoring point is selected at a certain interval, and the change rate of the rock mass strength between adjacent monitoring points is calculated. If the distance between two adjacent monitoring points is 2 meters and the difference in the rock mass strength is 12 MPa, then the change rate of the rock mass strength is 12 MPa÷2 m = 6 MPa / m, which exceeds the set threshold of 5 MPa / m.

[0037] When it is detected that the change rate of the rock mass strength in adjacent areas exceeds the set gradient threshold, reduce the cutting spacing. The original cutting spacing may be 5 meters. When the situation of exceeding the threshold occurs, the cutting spacing is reduced to 3 meters. In this way, in the area where the geological conditions change greatly, the model can be divided in more detail to improve the analytical accuracy of the geological attributes of the sub-model.

[0038] At the same time, when determining the boundary of each sub-model, the included angle between the boundary and the geological structure line is maintained within a set range through spatial vector calculation. The geological structure line can be determined according to the previously obtained geological radar detection data, core drilling data, and secondary supplementary exploration data. For example, by analyzing the strike and distribution of the rock mass joints, the general direction of the geological structure line is determined.

[0039] For the boundary of each sub-model, calculate the included angle between it and the geological structure line. The dot product and cross product operations of spatial vectors are used to calculate the included angle. The set range of the included angle between the boundary and the geological structure line is 30° - 60°. If the calculated included angle is not within this range, adjust the position of the boundary of the sub-model so that the included angle meets the requirements. In this way, acute-angle intersection between the cutting surface and the rock formation strike can be avoided, ensuring that the sub-model can more accurately reflect the geological structure.

[0040] Step S120: Perform gradient analysis on the longitudinal rock formation attributes and axial rock formation attributes of each of the sub-models to generate rock formation transition feature data of the sub-models, where the rock formation transition feature data includes the attribute change rate between adjacent rock formations and the spatial distribution of the transition region.

[0041] In this embodiment, after the cutting of the sub-models is completed, in-depth gradient analysis is performed on the longitudinal rock formation attributes and axial rock formation attributes of each sub-model to generate data reflecting the rock formation transition features. In this underground construction project, the geological conditions of the sub-models are complex and diverse, and the rock formation attributes in different sub-models are significantly different. Therefore, it is necessary to accurately analyze the rock formation transition features.

[0042] Step S121: Extract the sequence of rock mass quality indicators along the axis direction of the cavern. The sequence of rock mass quality indicators is generated by comprehensively calculating the borehole wave velocity test, the rock mass integrity coefficient, and the RQD value. Perform first-order difference calculation on the sequence of rock mass quality indicators to obtain the quality indicator change gradient between adjacent rock formations, and mark the section where the absolute value of the gradient of the quality indicator change gradient exceeds the set threshold as the mutation region.

[0043] To extract the sequence of rock mass quality indicators along the axis direction of the cavern, first perform a borehole wave velocity test. In each sub-model, a set number of boreholes are arranged along the axis direction of the cavern, and a sonic wave tester is placed in the boreholes. By emitting sonic waves into the boreholes, the propagation velocity of the sonic waves in the rock mass is measured. For example, in a sub-model, a borehole is arranged every 2 meters along the axis of the cavern, and a total of 10 boreholes are arranged, and the propagation velocity of the sonic waves at each borehole position is measured in turn.

[0044] At the same time, calculate the rock mass integrity coefficient. The rock mass integrity coefficient is determined by the ratio of the longitudinal wave velocity of the rock mass to the longitudinal wave velocity of the rock block. The longitudinal wave velocity of the rock block taken out from the borehole is tested in the laboratory, and combined with the longitudinal wave velocity of the rock mass obtained from the borehole wave velocity test, the rock mass integrity coefficient is calculated. For example, if the longitudinal wave velocity of the rock mass at a certain borehole position is 3000 m / s and the corresponding longitudinal wave velocity of the rock block is 4000 m / s, then the rock mass integrity coefficient at this position is 3000 m / s ÷ 4000 m / s = 0.75.

[0045] In addition, obtain the RQD value (rock quality designation). The RQD value is determined by statistically calculating the ratio of the cumulative length of the core segments with a length greater than 10 cm in the borehole to the footage of the borehole. During the core drilling process, the length of each core segment is carefully measured, and the cumulative length of the core segments with a length greater than 10 cm is statistically calculated. For example, if the footage of a certain borehole is 10 meters and the cumulative length of the core segments with a length greater than 10 cm is 7 meters, then the RQD value at this borehole position is 7 meters ÷ 10 meters = 70%.

[0046] Based on the comprehensive borehole wave velocity test, rock mass integrity coefficient, and RQD value, the rock mass quality index at each borehole location is calculated by the weighted average method. Assume that the weight of the borehole wave velocity test result is 0.4, the weight of the rock mass integrity coefficient is 0.3, and the weight of the RQD value is 0.3. The score of the borehole wave velocity test at a certain borehole location is converted to 80 points, the rock mass integrity coefficient is 0.75 (converted to 75 points), and the RQD value is 70 points. Then, the rock mass quality index at this location is 80×0.4 + 75×0.3 + 70×0.3 = 76.5 points.

[0047] Arrange the rock mass quality indexes at each borehole location along the axis direction of the cavern in sequence to form a rock mass quality index sequence. Perform a first-order difference calculation on this sequence, that is, calculate the difference between adjacent two rock mass quality indexes. For example, if the rock mass quality index sequence is [76.5, 78, 82, 79,...], then the first-order difference sequence is [78 - 76.5, 82 - 78, 79 - 82,...] = [1.5, 4, -3,...].

[0048] Set the threshold of the quality index change gradient to 3 points. Mark the section where the absolute value in the first-order difference sequence exceeds 3 points as the mutation area. For example, in the above first-order difference sequence, the absolute values of 4 and -3 both exceed 3 points, and the corresponding sections of the change in the rock mass quality index are marked as the mutation areas.

[0049] Step S122: Extract the distribution characteristics of the rock mass elastic modulus in the direction perpendicular to the axis of the cavern, convert the discrete borehole elastic modulus data in the distribution characteristics of the rock mass elastic modulus into a continuous spatial distribution matrix through the finite element interpolation algorithm, and calculate the Gaussian curvature of the spatial distribution matrix to generate a modulus change rate distribution map reflecting the spatial change rate of the rock mass stiffness.

[0050] In this embodiment, within each sub-model, extract the distribution characteristics of the rock mass elastic modulus in the direction perpendicular to the axis of the cavern. By conducting elastic modulus tests in boreholes at different positions within the sub-model, discrete borehole elastic modulus data is obtained. For example, in the sub-model, boreholes are arranged at intervals of 3 meters in the direction perpendicular to the axis of the cavern, and a total of 8 boreholes are arranged to measure the rock mass elastic modulus at each borehole location.

[0051] After obtaining the discrete borehole elastic modulus data, use the finite element interpolation algorithm to convert it into a continuous spatial distribution matrix. The finite element interpolation algorithm can estimate the elastic modulus of unmeasured positions based on the elastic modulus data at borehole positions and their spatial relationships. For example, for an unmeasured position located in the middle of multiple boreholes, the elastic modulus value of this position is estimated by weighted averaging according to the elastic modulus data of the surrounding boreholes in inverse proportion to the distance.

[0052] Divide the space within the sub-model into a set number of small grids, with each small grid corresponding to an element in the matrix. Fill the estimated elastic modulus values into the corresponding positions in the matrix to form a continuous spatial distribution matrix. The number of rows and columns of the spatial distribution matrix is determined according to the size of the sub-model and the density of the grid division.

[0053] Then, perform Gaussian curvature calculation on the generated spatial distribution matrix. Gaussian curvature is an index describing the degree of surface curvature, and in this scenario, it is used to reflect the spatial change rate of the rock mass stiffness. By calculating each element in the spatial distribution matrix and its adjacent elements, the Gaussian curvature value at each position is obtained.

[0054] For example, for an element in the spatial distribution matrix, consider the values of its four adjacent elements above, below, left, and right, and calculate the Gaussian curvature at this position through a common formula in the prior art. Combine the Gaussian curvature values at all positions to generate a modulus change rate distribution map reflecting the spatial change rate of the rock mass stiffness. In the modulus change rate distribution map, different colors or grayscales represent different Gaussian curvature values, thus visually showing the spatial change of the rock mass elastic modulus.

[0055] Step S123: Perform spatial overlay analysis on the position coordinates of the mutation region and the modulus change rate distribution map to determine the main gradient direction of the attribute change rate and the geometric shape of the transition region in the strata transition characteristic data. Among them, the geometric shape is extracted as a polygon patch through the boundary tracing algorithm, and its area, perimeter, and relative position relationship with the sub-model boundary are calculated.

[0056] Perform spatial overlay analysis on the position coordinates of the previously marked mutation region and the modulus change rate distribution map. The position coordinates of the mutation region can be determined by the position information of the drill holes and the coordinate system of the chamber axis. Thus, project the above coordinates onto the modulus change rate distribution map to find the corresponding position of the mutation region in the map.

[0057] By analyzing the distribution of the mutation region on the modulus change rate distribution map, determine the main gradient direction of the attribute change rate in the strata transition characteristic data. The main gradient direction represents the direction in which the rock mass attributes change the fastest. For example, observe the gradient change trend of the mutation region on the modulus change rate distribution map. If the increase or decrease of the modulus change rate is most obvious in a certain direction, then this direction is the main gradient direction.

[0058] When using the boundary tracking algorithm to extract the geometric shape of the transition region, the boundary tracking algorithm can start from the edge points of the mutation region in the modulus change rate distribution map and sequentially track adjacent points according to the set rules to gradually outline the boundary of the transition region. For example, starting from an edge point, check the points in its eight adjacent directions (up, down, left, right, upper left, lower left, upper right, lower right). If an adjacent point also belongs to the mutation region, add it to the boundary point sequence, and then use this point as the new starting point to continue tracking until returning to the starting point to form a closed boundary.

[0059] Connecting the above boundary points forms a polygon patch to represent the geometric shape of the transition region. Next, it is necessary to calculate the area of this polygon patch. Among them, the segmentation method can be used to divide the polygon into multiple triangles, calculate the area of each triangle and then sum them to obtain the area of the polygon. For example, for a quadrilateral polygon, one of its diagonals can be connected to divide it into two triangles. When calculating the area of a triangle, given the coordinates of the three vertices of the triangle, the method of vector cross product can be used to calculate. Suppose the coordinates of the three vertices of the triangle are A(x1, y1), B(x2, y2), and C(x3, y3) respectively, then the area of the triangle is 0.5 times the absolute value of the cross product of vector AB and vector AC. By sequentially calculating the area of each triangle after segmentation and summing them, the total area of the polygon patch can be obtained.

[0060] When calculating the perimeter of the polygon patch, just add the lengths of each side of the polygon. The length of each side can be calculated according to the distance formula between two points, that is, for two points P(x1, y1) and Q(x2, y2), the distance between them is the square root of the sum of the square of the difference in the abscissas of the two points and the square of the difference in the ordinates. Without using the formula editor, it is described as: first calculate the sum of the square of the difference in the abscissas of the two points and the square of the difference in the ordinates, and then take the square root of this sum. Sequentially calculate the lengths of each side of the polygon and accumulate them to obtain the perimeter.

[0061] When determining the relative position relationship between the polygon patch and the boundary of the sub-model, it is necessary to compare the vertex coordinates of the polygon with the coordinate range of the sub-model boundary. For example, the coordinate range of the sub-model on the plane is x from xmin to xmax, and y from ymin to ymax. For each vertex of the polygon, judge whether its x coordinate and y coordinate are within the coordinate range of the sub-model. If all vertices of the polygon are within the sub-model boundary, it means that the transition region is completely inside the sub-model; if some vertices are on the sub-model boundary, the transition region intersects with the sub-model boundary; if there are vertices outside the sub-model boundary, the transition region partially exceeds the sub-model range. Parameters such as the overlapping area or the exceeding area between the transition region and the sub-model boundary can also be further calculated to more accurately describe their relative position relationship.

[0062] Step S130: According to the rock formation transition feature data, dynamically identify the surrounding rock types of each sub-model, map and match the surrounding rock identification results with the spatial coordinates of the preset protected object to obtain a mapping and matching relationship, and collect the geological parameter set of the sub-model in real time based on the mapping and matching relationship.

[0063] Step S131: Dynamically identify the surrounding rock types of each sub-model according to the rock formation transition feature data, including: In this embodiment, the rock formation transition feature data can be input into a pre-trained surrounding rock classification model to output a surrounding rock category probability distribution vector. Among them, the surrounding rock classification model is constructed based on a convolutional neural network architecture, and its training data includes sample rock formation transition feature data of different surrounding rock types in historical projects.

[0064] In this underground construction project, the rock formation transition feature data of each sub-model has been obtained, including multi-dimensional information such as the main gradient direction of the attribute change rate and the geometric shape of the transition area. Organize the above rock formation transition feature data into a format suitable for input into the model. For example, represent the main gradient direction in vector form and splice the geometric shape features (area, perimeter, etc.) of the transition area as numerical features.

[0065] The pre-trained surrounding rock classification model is constructed based on a convolutional neural network architecture. This surrounding rock classification model uses a large amount of sample rock formation transition feature data of different surrounding rock types in historical projects during the training stage. During the training process, the surrounding rock classification model continuously adjusts its own weight parameters to learn the rock formation transition feature patterns corresponding to different surrounding rock types.

[0066] When the rock formation transition feature data of the current sub-model is input into the pre-trained surrounding rock classification model, the surrounding rock classification model performs a series of operations such as convolution, pooling, and fully connected. The convolution operation extracts local features in the rock formation transition feature data, the pooling operation reduces the dimension of the local features, and the fully connected layer integrates the extracted local features and outputs the result. Finally, a surrounding rock category probability distribution vector is output. Each element of this surrounding rock category probability distribution vector represents the probability that the surrounding rock of this sub-model belongs to a certain set type. For example, the vector [0.1, 0.2, 0.7] means that the probability that the surrounding rock of this sub-model belongs to the first type is 0.1, the probability that it belongs to the second type is 0.2, and the probability that it belongs to the third type is 0.7.

[0067] Step S132: Generate the surrounding rock type identification result of the surrounding rock type label according to the probability distribution vector of the surrounding rock category, and generate the corresponding three-dimensional surrounding rock type distribution heat map in combination with the spatial position coordinates of the sub-model. The rendering rule of the three-dimensional surrounding rock type distribution heat map is: assign a set color code to different surrounding rock categories, and adjust the transparency according to the probability value of the implemented surrounding rock type label to reflect the classification confidence level.

[0068] In this embodiment, according to the output probability distribution vector of the surrounding rock category, the surrounding rock type corresponding to the element with the largest probability value can be selected as the surrounding rock type label of the sub-model. For example, in the above probability distribution vector [0.1, 0.2, 0.7], the surrounding rock type corresponding to the maximum probability value 0.7 is the surrounding rock type label of the sub-model, and thus the surrounding rock identification result of each sub-model is obtained.

[0069] Then, in combination with the spatial position coordinates of the sub-model, generate the corresponding three-dimensional surrounding rock type distribution heat map. In three-dimensional space, each sub-model has its specific coordinate position. Assign a set color code to different surrounding rock categories. For example, set the hard granite surrounding rock to red and the fractured sandstone surrounding rock to blue, etc. At the same time, adjust the transparency of the sub-model in the heat map according to the probability value of the surrounding rock type label. The higher the probability value, the lower the transparency, indicating a higher classification confidence level; the lower the probability value, the higher the transparency, indicating a lower classification confidence level.

[0070] For example, the probability value corresponding to the surrounding rock type label of a sub-model is 0.9, and the sub-model is displayed in the heat map with a lower transparency for the corresponding color; the probability value of another sub-model is 0.3, and it is displayed in the heat map with a higher transparency for the color, so that the confidence level of the surrounding rock classification of each sub-model can be intuitively reflected.

[0071] Step S133: Optimize the boundary of the area where the difference between the surrounding rock type labels of adjacent sub-models in the surrounding rock type distribution heat map exceeds a set value. Among them, the method of boundary optimization includes: extracting the rock layer transition feature data of the boundary area for secondary classification. If the secondary classification result is consistent with the surrounding rock identification result, the original boundary is maintained; otherwise, the sub-model cutting line is adjusted based on the secondary classification result to generate the corrected surrounding rock identification result.

[0072] In this embodiment, in the generated three-dimensional surrounding rock type distribution heat map, check the surrounding rock type labels of adjacent sub-models. Among them, a difference threshold can be set. For example, when the surrounding rock type labels of adjacent sub-models are completely different, it is considered that the difference exceeds the set value. For the above area where the difference exceeds the set value, determine it as the boundary area.

[0073] Then, extract the rock stratum transition feature data of the boundary area and input it again into the pre-trained surrounding rock classification model for secondary classification. For example, for a sub-model in the boundary area, after organizing its rock stratum transition feature data, re-enter it into the surrounding rock classification model to obtain a new probability distribution vector of the surrounding rock category and determine a new surrounding rock category label.

[0074] Next, compare the secondary classification result with the previous surrounding rock identification result. If the secondary classification result is consistent with the original surrounding rock identification result, it indicates that the original classification is accurate and the original boundary remains unchanged. If the secondary classification result is different from the original surrounding rock identification result, then adjust the sub-model cutting line based on the secondary classification result. For example, originally the surrounding rock category labels of two adjacent sub-models are different, but after secondary classification, it is found that they should belong to the same type, then adjust the cutting line between these two sub-models so that they are merged into a sub-model with the same surrounding rock type. Through such a boundary optimization process, a corrected surrounding rock identification result is generated.

[0075] Step S134: Map and match the surrounding rock identification result with the spatial coordinates of the preset protected object to obtain a mapping and matching relationship, and based on the mapping and matching relationship, collect the set of geological parameters of the sub-model in real time, including: In this embodiment, determine the corresponding three-dimensional influence range according to the geographical coordinates of the preset protected object. The three-dimensional influence range is calculated and generated by a geological disturbance propagation model. The inputs of the geological disturbance propagation model include the structural type of the protected object, the foundation burial depth, and the rock mass stress distribution parameters.

[0076] In an underground construction project, the preset protected object may include surrounding buildings, underground pipelines, etc. First, obtain the geographical coordinates of the above-mentioned preset protected object, which can be obtained by means such as a geographic information system (GIS).

[0077] The geological disturbance propagation model is used to calculate the three-dimensional influence range of the protected object. The inputs of the geological disturbance propagation model include the structural type of the protected object (such as frame structure, brick-concrete structure, etc.), the foundation burial depth (for example, the foundation burial depth of a certain building is 5 meters), and the rock mass stress distribution parameters. The rock mass stress distribution parameters can be obtained from previous geological exploration and monitoring data, such as the rock mass stress values measured at different positions.

[0078] The geological disturbance propagation model can consider various factors to calculate the three-dimensional influence range. For protection objects of different structural types, their abilities to resist geological disturbances are different. For example, frame structures are relatively more capable of withstanding certain deformations. The foundation burial depth also affects the propagation of geological disturbances. The influence range on a foundation with a deeper burial depth may be relatively smaller. The geological disturbance propagation model can, based on the above input parameters, determine the three-dimensional spatial range around the protection object that may be affected by geological disturbances through calculation and simulation. For example, for a frame structure building with a foundation burial depth of 5 meters, the geological disturbance propagation model calculates that under the influence of geological disturbances that may be generated during the construction of an underground cavity, the three-dimensional influence range is a cylindrical space centered on the building, with a radius of 20 meters and a height from the ground to 10 meters underground.

[0079] Step S135: Project the three-dimensional influence range into the three-dimensional model, determine the set of sub-models that intersect with it through spatial topological relationship analysis, and establish the mapping and matching relationship between the preset protection object and the sub-models.

[0080] In this embodiment, the calculated three-dimensional influence range of the protection object is projected into the three-dimensional model of the underground cavity established previously. During the projection process, the spatial coordinates of the three-dimensional influence range can be aligned with the coordinate system of the three-dimensional model to determine its specific position in the three-dimensional model.

[0081] Through spatial topological relationship analysis, check whether the three-dimensional influence range intersects with the sub-models in the three-dimensional model. Spatial topological relationship analysis will compare the boundaries of the three-dimensional influence range with the boundaries of the sub-models to determine whether there is an overlapping part. For example, for a sub-model and the three-dimensional influence range, if the boundary of the three-dimensional influence range intersects with the boundary of the sub-model, or the three-dimensional influence range completely contains the sub-model, or the sub-model completely contains the three-dimensional influence range, it is considered that they intersect.

[0082] Form a set of all sub-models that intersect with the three-dimensional influence range, and establish the mapping and matching relationship between the preset protection object and the above sub-models. For example, for a protection object, if the determined sub-models that intersect with it are Sub-model A, Sub-model B, and Sub-model C, then establish the mapping relationship between the protection object and Sub-models A, B, and C, indicating that the geological conditions of the above sub-models may affect this protection object.

[0083] Step S136: Based on the mapping and matching relationship, collect in real time the set of geological parameters of the sub-models, and extract the surrounding rock identification results of the set of sub-models, and calculate the geological risk index of the area where the preset protection object is located. The specific calculation logic is as follows: Based on the stability coefficient corresponding to the surrounding rock category label, the groundwater pressure value and the in-situ stress concentration degree in the set of geological parameters of the sub-models, construct a weighted evaluation model and output the risk index level.

[0084] In this embodiment, according to the established mapping matching relationship, the geological parameter set related to the preset protected object is collected in real time, which may include information such as groundwater pressure value, in-situ stress concentration degree, etc. The groundwater pressure value can be obtained by real-time monitoring with piezometers installed in the sub-model area. For example, the piezometer in a certain sub-model area shows that the groundwater pressure is 0.2 MPa. The in-situ stress concentration degree can be measured by stress sensors. For example, the in-situ stress concentration degree in a certain sub-model area is 1.5.

[0085] At the same time, the surrounding rock identification results of the above-mentioned sub-models are extracted to obtain the surrounding rock category labels of each sub-model. Different surrounding rock category labels correspond to different stability coefficients. For example, the stability coefficient of hard granite surrounding rock is 0.9, and the stability coefficient of fractured sandstone surrounding rock is 0.3.

[0086] A weighted evaluation model is constructed to calculate the geological risk index of the area where the preset protected object is located. Assume that the weight of the stability coefficient is 0.4, the weight of the groundwater pressure value is 0.3, and the weight of the in-situ stress concentration degree is 0.3. For a set of sub-models related to the protected object, first calculate the risk score of each sub-model. For example, the stability coefficient corresponding to the surrounding rock category label of a sub-model is 0.7, the groundwater pressure value is 0.15 MPa, and the in-situ stress concentration degree is 1.2. Then the risk score of this sub-model = 0.7×0.4 + 0.15×0.3 + 1.2×0.3 = 0.625.

[0087] The risk scores of all sub-models related to the protected object are weighted and averaged (the weights here can be determined according to factors such as the spatial distance between the sub-model and the protected object) to obtain the geological risk index of the area where the preset protected object is located. According to the magnitude of the risk index, different risk index levels are divided. For example, the risk index is in the range of 0 - 0.3 for the low-risk level, 0.3 - 0.6 for the medium-risk level, and 0.6 - 1 for the high-risk level.

[0088] Step S 137: When the risk index level exceeds the preset threshold, a dynamic monitoring link is established between the preset protected object and the corresponding sub-model. The dynamic monitoring link collects the geological parameter change data in real time through the deployment of fiber optic strain sensors, piezometers, and microseismic monitoring equipment, and transmits the geological parameter change data to the central processing system for early warning analysis.

[0089] Among them, a preset threshold for the risk index level can be set. For example, the upper limit of the medium-risk level is 0.6. When the calculated geological risk index level of the area where the preset protected object is located exceeds this threshold, it indicates that the protected object faces a relatively high geological risk and needs to strengthen monitoring.

[0090] Then, a dynamic monitoring link is established between the pre-defined protected object and the corresponding sub-model. Fiber optic strain sensors, piezometers, and microseismic monitoring equipment are deployed in the corresponding sub-model areas. Fiber optic strain sensors can monitor rock strain in real time. For example, when the rock undergoes slight deformation, the optical signal changes detected by the fiber optic strain sensor. Analysis of the optical signal reveals the rock strain value. The piezometers continue to monitor groundwater pressure changes in real time. Microseismic monitoring equipment can detect subtle vibrations within the rock mass, such as those caused by crack expansion. These signals can be captured by the microseismic monitoring equipment.

[0091] The monitoring equipment collects real-time data on geological parameter changes, such as strain, groundwater pressure changes, and microseismic signals, and transmits it to a central processing system via wired or wireless communication. The central processing system analyzes this data and compares it with pre-set warning thresholds. For example, if strain exceeds a pre-set threshold, if groundwater pressure changes exceed a certain range, or if the frequency and intensity of microseismic signals reach a certain level, the central processing system triggers a warning signal, alerting personnel to take appropriate measures.

[0092] Step S140: triggering a safety warning signal according to the numerical range of the geological parameter set, and generating a set of emergency response plans associated with the safety warning signal.

[0093] Step S141: normalizing the actual measured value of each geological parameter in the geological parameter set with its corresponding safety threshold to generate a deviation coefficient for each parameter.

[0094] In this embodiment, in underground construction projects, the geological parameter set includes multiple parameters, such as groundwater pressure, geostress, and rock strain. Each geological parameter has a corresponding safety threshold, which is predetermined based on geological conditions, engineering design requirements, and previous experience. For example, the safety threshold for groundwater pressure is set at 0.3 MPa, the safety threshold for geostress is set at 2.0 MPa, and the safety threshold for rock strain is set at 0.001.

[0095] Therefore, for each geological parameter in the geological parameter set, its actual measured value is normalized against the corresponding safety threshold. The purpose of normalization is to unify the actual measured values of different parameters onto a comparable scale. Taking groundwater pressure as an example, assuming that the groundwater pressure value obtained in a certain actual measurement is 0.35 MPa, during normalization, the difference between the actual measured value and the safety threshold is first calculated: 0.35 MPa - 0.3 MPa = 0.05 MPa. This difference is then divided by the safety threshold to obtain the deviation coefficient: 0.05 MPa ÷ 0.3 MPa ≈ 0.167.

[0096] The same method is used to normalize other parameters in the geological parameter set, such as in-situ stress and rock strain, to obtain a deviation coefficient for each parameter. For example, if the in-situ stress value measured is 2.2 MPa, its deviation coefficient is (2.2 MPa - 2.0 MPa) ÷ 2.0 MPa = 0.1; if the rock strain value measured is 0.0012, its deviation coefficient is (0.0012 - 0.001) ÷ 0.001 = 0.2.

[0097] Step S142: performing weighted calculation on the deviation coefficient based on a preset deviation weight matrix to generate a comprehensive risk value.

[0098] The preset deviation weight matrix is determined based on the importance of each geological parameter to underground cavern safety. For example, for underground caverns, changes in geostress and groundwater pressure have a greater impact on safety, while rock strain has a relatively smaller impact. Therefore, the deviation weight for geostress can be set to 0.4, the deviation weight for groundwater pressure to 0.3, and the deviation weight for rock strain to 0.3.

[0099] Therefore, the deviation coefficient of each geological parameter is multiplied by the corresponding deviation weight, and then these products are added together to obtain the comprehensive risk value. Taking the deviation coefficient calculated above as an example, the comprehensive risk value = 0.1 × 0.4 + 0.167 × 0.3 + 0.2 × 0.3 = 0.1501.

[0100] Step S143: When the comprehensive risk value exceeds a first threshold, a first warning signal is triggered and a preliminary treatment plan set is generated. The preliminary treatment plan set includes a grouting reinforcement range, support parameter adjustment instructions, and monitoring frequency increase parameters.

[0101] The first threshold value can be set to 0.1. When the calculated comprehensive risk value of 0.1501 exceeds the first threshold value, a first warning signal is triggered. The first warning signal can be issued by an audible and visual alarm device to alert relevant personnel to the increase in geological risk.

[0102] At the same time, a preliminary set of treatment options is generated. To determine the scope of grouting reinforcement, areas of potential safety hazards are analyzed based on the changes in geological parameters and the distribution of submodels. For example, if the deviation coefficient of ground stress and groundwater pressure in a submodel area is large, then grouting reinforcement may be required in the area containing that submodel and its adjacent submodels. After determining the scope of grouting reinforcement, a specific grouting plan is developed, including parameters such as grouting material selection, grouting pressure, and grouting volume.

[0103] The support parameter adjustment instruction can adjust the support structure of the underground cavern according to the geological risk situation. For example, if the in-situ stress increases, it may be necessary to increase the length and density of the rock bolts or improve the strength grade of the support material. According to the specific changes in the comprehensive risk value and geological parameters, the appropriate support parameter adjustment value is calculated.

[0104] The monitoring frequency increase parameter is for closer monitoring of geological parameters. When determining the monitoring frequency increase parameter, it is adjusted based on the deviation coefficient of geological parameters and the comprehensive risk value. If the deviation coefficients of in-situ stress and groundwater pressure are relatively high, it indicates that the changes in these two parameters are more dangerous and the monitoring frequency for them needs to be increased. For example, originally the in-situ stress and groundwater pressure were monitored once every 2 hours. When the comprehensive risk value exceeds the first threshold, the monitoring frequency is increased to once every 1 hour. For rock mass strain, if its deviation coefficient is relatively small, the monitoring frequency can be increased from once every 4 hours to once every 2 hours.

[0105] Step S144: When the comprehensive risk value exceeds the second threshold, trigger the second warning signal and activate the emergency disposal process. The emergency disposal process includes a construction stop instruction, a drainage system start instruction, and a multi-objective optimization plan generation instruction.

[0106] Among them, the second threshold can be set to 0.2. If the comprehensive risk value continues to rise and exceeds this second threshold of 0.2, the situation is more critical at this time, and the second warning signal is triggered. The second warning signal is stronger than the first warning signal and may notify all relevant persons in charge and construction personnel of the project through various means such as high-decibel alarm sounds, flashing red lights, text messages, and emails.

[0107] Immediately activate the emergency disposal process. First, issue a construction stop instruction, and quickly convey the information to stop all construction activities through devices such as the broadcast system and walkie-talkies at the construction site. The construction personnel must immediately stop the operations in their hands and evacuate to the safe area in an orderly manner. This is to avoid more serious safety accidents caused by continuing construction under unstable geological conditions.

[0108] At the same time, start the drainage system. The start of the drainage system needs to be operated according to the specific layout of the underground cavern and the groundwater situation. If there are multiple drainage zones in the underground cavern, it is necessary to first determine which zone has too high groundwater pressure, and then give priority to starting the drainage pump in that zone. For example, based on the data fed back by the water level sensors and pressure sensors installed in different zones, it is judged that the groundwater pressure in a certain zone has reached a dangerous value, and the drainage pump in that zone is immediately turned on for drainage operations to reduce the groundwater pressure and minimize the impact on the cavern structure.

[0109] In addition, generate multi-objective optimization scheme generation instructions. This requires comprehensive consideration of multiple factors, such as construction progress, cost, safety, etc. By collecting current geological parameters, construction progress data, and previous construction plans and other information, optimization algorithms are used to generate multi-objective optimization schemes. For example, considering both ensuring construction safety and minimizing the impact on construction progress, the algorithm will optimize aspects such as adjusting the support scheme and changing the construction sequence. If the geological conditions deteriorate, the original plan of excavating a large area at one time may be adjusted to excavate in multiple small areas, while increasing temporary support measures to ensure construction safety and minimize the delay of the construction period as much as possible.

[0110] Step S145: According to the first warning signal or the second warning signal, screen the disposal cases matching the current geological parameter set from the historical database, and extract the schemes meeting the preset safety indicators in the disposal cases as the candidate scheme set.

[0111] In this embodiment, after the first warning signal or the second warning signal is triggered, the disposal cases matching the current geological parameter set can be automatically screened from the historical database. The historical database stores the disposal cases of previous similar projects under different geological conditions, including information such as geological parameter conditions, adopted disposal measures, and disposal effects.

[0112] Then, compare the current geological parameter set, such as the parameter values of groundwater pressure, ground stress, rock mass strain, etc., with the cases in the historical database. During the comparison process, a parameter matching range will be set. For example, the matching range of groundwater pressure is ±10% of the current value, and the matching range of ground stress is ±15%, etc. For each historical case, check whether its geological parameters are within the matching range of the current parameters.

[0113] For example, the current groundwater pressure is 0.35 MPa, the ground stress is 2.2 MPa, and the rock mass strain is 0.0012. Searching in the historical database, a case is found where the groundwater pressure is 0.33 MPa, the ground stress is 2.3 MPa, and the rock mass strain is 0.0011. Since the above parameters are all within the matching range of the current parameters, this case is screened out.

[0114] After screening out the matching cases, extract the schemes meeting the preset safety indicators. The preset safety indicators include requirements for the stability of the support structure, requirements for groundwater control, etc. For example, it is required that the support structure can withstand a certain ground stress after disposal, and the groundwater pressure should be reduced to a safe range, etc. Only the schemes meeting the above safety indicators will be used as candidate schemes.

[0115] Step S146: Verify the construction conditions of the candidate solution set, and eliminate infeasible solutions according to the available status of the current construction equipment, the material inventory, and the construction period constraint conditions.

[0116] In this embodiment, after obtaining the candidate solution set, it is necessary to verify the construction conditions of the above candidate solutions. Specifically, first consider the available status of the current construction equipment. Different disposal solutions may require different construction equipment. For example, the grouting reinforcement solution requires equipment such as grouting pumps, and the solution of adding a support structure requires equipment such as anchor drills. Check whether the above-mentioned equipment existing at the construction site is in an available state and whether there is a sufficient quantity to meet the implementation requirements of the solution.

[0117] For example, if a certain candidate solution requires 3 grouting pumps to operate simultaneously, but there are only 2 available grouting pumps at the construction site currently, and the other one is under repair and cannot be put into use in the short term, then this solution does not meet the requirements of the available status of the construction equipment and will be eliminated.

[0118] Then check the material inventory. The implementation of the disposal solution often requires the consumption of certain materials, such as grouting materials, anchor bolts, steel, etc. Check whether the current inventory quantity of materials can meet the requirements of the solution. If a certain solution requires 100 tons of grouting materials, but the current inventory is only 50 tons, and it is impossible to replenish enough materials in the short term, then this solution will also be eliminated.

[0119] It is also necessary to consider the construction period constraint conditions. Each project has certain construction period requirements. If the implementation of a certain candidate solution will cause the construction period to be delayed too long and exceed the allowable range of the project, it will also be considered an infeasible solution. For example, if the project requires a certain stage of construction to be completed within 1 month, and a certain candidate solution requires 2 months to complete, then this solution does not meet the construction period constraint conditions and will be eliminated.

[0120] Step S147: Sort the verified candidate solution set according to the implementation efficiency, cost consumption, and risk suppression effect to generate the final emergency disposal solution set.

[0121] In this embodiment, after the construction condition matching verification, a feasible candidate solution set is obtained. Next, sort the above-mentioned feasible candidate solution set in terms of implementation efficiency, cost consumption, and risk suppression effect.

[0122] Specifically, in terms of implementation efficiency, evaluate the time required for each solution to achieve the expected effect from the start of implementation. For example, if Solution A is expected to be completed in 10 days and Solution B is expected to be completed in 15 days, then the implementation efficiency of Solution A is relatively high. The implementation efficiency can be determined by analyzing factors such as the time arrangement of each construction step in the solution, the allocation of equipment and personnel, etc.

[0123] In terms of cost consumption, calculate the expenses required during the implementation of each plan, including equipment rental costs, material procurement costs, labor costs, etc. For example, if the total cost of Plan A is 500,000 yuan and the total cost of Plan B is 800,000 yuan, then the cost consumption of Plan A is relatively low.

[0124] In terms of risk suppression effect, evaluate the effect of each plan on reducing geological risks through simulation analysis or referring to historical cases. For example, after the implementation of Plan A, it is expected that the comprehensive risk value can be reduced to below 0.1, and after the implementation of Plan B, it is expected that the comprehensive risk value can be reduced to 0.12. Then the risk suppression effect of Plan A is better.

[0125] Then, set weights for the above three aspects respectively. For example, the weight of implementation efficiency is 0.3, the weight of cost consumption is 0.3, and the weight of risk suppression effect is 0.4. Score the implementation efficiency, cost consumption, and risk suppression effect of each plan, and then sort them according to the weighted scores. The plan with the highest score is ranked first, and the plan with the lowest score is ranked last.

[0126] Finally, use the sorted set of candidate plans as the final set of emergency response plans. In practical applications, the plan ranked first can be preferentially selected for implementation. If problems are encountered during the implementation of this plan, then consider the subsequent plans in turn.

[0127] Through the above steps, the embodiment of the present invention realizes the efficient and accurate analysis of potential safety hazards in underground caverns, significantly improving the safety and reliability of underground cavern projects. Specifically, first, based on the three-dimensional model data of the underground cavern, combined with the dynamically adjusted geological condition segmentation rules, the complex three-dimensional geological model is finely cut into multiple sub-models, which not only considers the gradient change characteristics of the rock strata but also ensures that each sub-model can accurately reflect the subtle differences in local geological conditions. On this basis, further gradient analysis is carried out on the longitudinal rock stratum attributes and axial rock stratum attributes of each sub-model to generate rock stratum transition characteristic data including the attribute change rate between adjacent rock strata and the spatial distribution of the transition area, which not only reveals the internal connection and change law between rock strata but also provides a key basis for the dynamic identification of surrounding rock types. By mapping and matching the surrounding rock identification results with the spatial coordinates of the preset protected object, a set of geological parameters closely related to the protected object can be collected in real time, thus realizing the accurate positioning and real-time monitoring of potential safety hazards. Finally, according to the numerical interval of the set of geological parameters, the method can intelligently trigger safety warning signals and generate an associated set of emergency response plans, improving the response speed to potential safety hazards. In summary, a comprehensive, in-depth analysis and efficient response to potential safety hazards in underground caverns are realized, which helps the safe construction and operation of underground cavern projects.

[0128] Step S210: The method further includes a dynamic update step for the sub-model: Step S211: Real-time collect the vibration spectrum data and hydraulic propulsion pressure data of the tunneling equipment, and extract the main vibration frequency characteristics and pressure fluctuation amplitude.

[0129] During the tunneling process of the underground chamber, relevant data of the tunneling equipment can be collected in real time. For the collection of vibration spectrum data, vibration sensors are installed at key parts of the tunneling equipment, such as the cutter head and the fuselage. The above vibration sensors can monitor the vibration of the equipment during operation in real time and convert the vibration signal into an electrical signal and transmit it to the data acquisition system.

[0130] The data acquisition system can process the vibration signal and convert the vibration signal in the time domain into vibration spectrum data in the frequency domain through methods such as Fourier transform. For example, analyze the vibration signal within a period of time to obtain the vibration amplitudes of different frequency components. Extract the main vibration frequency characteristics from the above vibration spectrum data, that is, the frequency component with the largest vibration amplitude. For example, after analysis, it is found that the component with a frequency of 50Hz in the vibration spectrum has the largest vibration amplitude, then 50Hz is the main vibration frequency at this moment.

[0131] At the same time, collect the hydraulic propulsion pressure data. Install a pressure sensor in the hydraulic system to measure the value of the hydraulic propulsion pressure in real time. The pressure sensor converts the pressure signal into an electrical signal and transmits it to the data acquisition system. The data acquisition system will record the hydraulic propulsion pressure values at different times and calculate the pressure fluctuation amplitude. The pressure fluctuation amplitude refers to the difference between the maximum value and the minimum value of the hydraulic propulsion pressure within a certain period of time. For example, within 1 minute, the maximum value of the hydraulic propulsion pressure is 20MPa, and the minimum value is 18MPa, then the pressure fluctuation amplitude is 20MPa - 18MPa = 2MPa.

[0132] Step S212: Input the main vibration frequency characteristics and pressure fluctuation amplitude into a preset rock mass strength-geological parameter mapping model to invert the predicted value of the rock mass strength and the fracture development density of the current excavation face.

[0133] Step S2121: Perform normalization processing on the main vibration frequency characteristics and pressure fluctuation amplitude respectively to generate a normalized main vibration frequency characteristic and a normalized pressure fluctuation amplitude, where the scaling factor of the normalization processing is dynamically adjusted according to the range difference of the historical vibration spectrum data and the hydraulic propulsion pressure.

[0134] In this embodiment, the purpose of performing normalization processing on the collected main vibration frequency characteristics and pressure fluctuation amplitude is to unify these two data with different ranges to a suitable scale for input into the preset rock mass strength-geological parameter mapping model.

[0135] The scaling factor for normalization is dynamically adjusted according to the range of historical vibration spectrum data and hydraulic propulsion pressure. The historical vibration spectrum data contains the vibration main frequency values at each moment during previous tunneling. By finding the maximum and minimum values among the above values, the range is calculated. For example, if the maximum value of the historical vibration main frequency is 80 Hz and the minimum value is 20 Hz, then the range is 80 Hz - 20 Hz = 60 Hz. Similarly, analyze the historical hydraulic propulsion pressure data, find the maximum and minimum values, and calculate the range.

[0136] For the current vibration main frequency feature, assuming its value is 60 Hz, during normalization, first calculate the difference between it and the historical minimum value, that is, 60 Hz - 20 Hz = 40 Hz, and then divide this difference by the range to obtain the normalized vibration main frequency feature as 40 Hz ÷ 60 Hz ≈ 0.67.

[0137] For the pressure fluctuation amplitude, assuming the current pressure fluctuation amplitude is 3 MPa, the maximum value of the historical hydraulic propulsion pressure is 25 MPa, and the minimum value is 15 MPa, and the range is 25 MPa - 15 MPa = 10 MPa. Then the normalized pressure fluctuation amplitude is (3 MPa) ÷ (10 MPa) = 0.3.

[0138] Step S2122: Align and splice the normalized vibration main frequency feature and the normalized pressure fluctuation amplitude according to the time stamp into a multi-dimensional feature vector, input it into the rock mass strength - geological parameter mapping model for forward calculation, and output the predicted value of the rock mass strength and the fracture development density of the current excavation face.

[0139] Step S21221: Extract the distribution weight of the vibration main frequency band and the time series change slope of the pressure fluctuation amplitude in the multi-dimensional feature vector to generate a frequency domain - time domain correlation feature matrix.

[0140] In this embodiment, the normalized vibration main frequency feature and the normalized pressure fluctuation amplitude can be aligned according to the time stamp, that is, to ensure that these two features at the same moment correspond. For example, at a certain moment t, the normalized vibration main frequency feature is 0.67 and the normalized pressure fluctuation amplitude is 0.3. Combine them to form a two-dimensional vector [0.67, 0.3]. As time goes by, a series of such two-dimensional vectors will be obtained, and they are spliced together to form a multi-dimensional feature vector.

[0141] Then, extract the distribution weights of the main vibration frequency bands from the multi-dimensional feature vector. Among them, the frequency band where the main vibration frequency is located can be divided into several sub-frequency bands, the number of occurrences of the main vibration frequency in each sub-frequency band is counted, and then the weight of each sub-frequency band is calculated. For example, the main vibration frequency band is divided into three sub-frequency bands of 20 - 40 Hz, 40 - 60 Hz, and 60 - 80 Hz. In a period of time, the main vibration frequency appears 10 times in the 20 - 40 Hz sub-frequency band, 20 times in the 40 - 60 Hz sub-frequency band, and 5 times in the 60 - 80 Hz sub-frequency band. Then the weights of these three sub-frequency bands are 10÷(10 + 20 + 5)≈0.29, 20÷(10 + 20 + 5)≈0.57, and 5÷(10 + 20 + 5)≈0.14 respectively.

[0142] Meanwhile, calculate the temporal change slope of the pressure fluctuation amplitude. By performing linear fitting on the normalized pressure fluctuation amplitudes at different times, the slope of the fitting line is obtained. For example, the normalized pressure fluctuation amplitudes recorded at 5 times are 0.2, 0.25, 0.3, 0.35, and 0.4 respectively, and the slope of the fitting line obtained through linear regression analysis is 0.05.

[0143] Combine the distribution weights of the main vibration frequency bands and the temporal change slope of the pressure fluctuation amplitude to generate a frequency-domain - time-domain correlation feature matrix. For example, combine the weights of the above three sub-frequency bands and the temporal change slope of the pressure fluctuation amplitude into a 4D vector [0.29, 0.57, 0.14, 0.05] as a row of the frequency-domain - time-domain correlation feature matrix. As time goes by, multiple such vectors will be obtained to form a matrix.

[0144] Step S21222: Input the frequency-domain - time-domain correlation feature matrix into the hidden units of the rock mass strength - geological parameter mapping model, and calculate the output vector of the hidden layer through the activation function.

[0145] In this embodiment, the generated frequency-domain - time-domain correlation feature matrix is input into the hidden units of the rock mass strength - geological parameter mapping model. The rock mass strength - geological parameter mapping model is a neural network model, and the hidden units are the intermediate layers in the model, which contain multiple neurons.

[0146] Each neuron will perform weighted summation on the input frequency-domain - time-domain correlation feature matrix, and then perform non-linear transformation through the activation function. The activation function can adopt functions such as the Sigmoid function and the ReLU function. Taking the Sigmoid function as an example, its role is to map the weighted summation result of the neuron to the range of 0 - 1.

[0147] For example, for a neuron in the hidden unit, it receives a row vector [0.29, 0.57, 0.14, 0.05] of the frequency-domain time-domain correlation feature matrix and assigns a weight to each element. Suppose the weights are 0.2, 0.3, 0.4, and 0.1 respectively. First, perform weighted summation, that is, 0.29×0.2 + 0.57×0.3 + 0.14×0.4 + 0.05×0.1 = 0.282. Then input this result into the Sigmoid function. The expression of the Sigmoid function is 1÷(1 + e to the negative x power), where x = 0.282, and calculate the output value of this neuron.

[0148] Perform such calculations for each neuron in the hidden unit to obtain the output vector of the hidden layer. For example, if there are 5 neurons in the hidden unit, the output vector of the hidden layer obtained after calculation is [0.4, 0.6, 0.3, 0.7, 0.2].

[0149] Step S21223: Input the output vector of the hidden layer into the rock mass strength regression branch and the fracture density regression branch respectively. The weight parameters of the rock mass strength regression branch are supervised and trained through historical borehole core strength data, and the weight parameters of the fracture density regression branch are supervised and trained through historical borehole fracture imaging data, and output the predicted value of the rock mass strength and the fracture development density of the current excavation face.

[0150] Input the output vector of the hidden layer into the rock mass strength regression branch and the fracture density regression branch respectively. The rock mass strength regression branch and the fracture density regression branch are two different branches of the rock mass strength-geological parameter mapping model, which are used to predict the rock mass strength and the fracture development density respectively.

[0151] The weight parameters of the rock mass strength regression branch are obtained through supervised training with historical borehole core strength data. In the historical data, there are strength test results of a large number of borehole core samples. Using the above historical borehole core strength data as the training target, adjust the weight parameters of the rock mass strength regression branch so that the model can predict the rock mass strength as accurately as possible according to the input output vector of the hidden layer.

[0152] The weight parameters of the fracture density regression branch are obtained through supervised training with historical borehole fracture imaging data. The historical borehole fracture imaging data records the fracture development conditions at different borehole positions. Using the above data as the training target, adjust the weight parameters of the fracture density regression branch so that the model can predict the fracture development density according to the input output vector of the hidden layer.

[0153] For example, input the hidden layer output vector [0.4, 0.6, 0.3, 0.7, 0.2] into the rock mass strength regression branch. After a series of linear transformations and calculations, output the predicted value of the rock mass strength of the current excavation face, which is assumed to be 120 MPa. Input this hidden layer output vector into the fracture density regression branch. After corresponding calculations, output the fracture development density of the current excavation face, which is assumed to be 5 fractures per cubic meter.

[0154] Step S213: Calculate the deviation rate between the predicted value of the rock mass strength and the designed value of the rock mass strength stored in the sub-model. At the same time, compare the pressure fluctuation amplitude with the preset pressure threshold. If the deviation rate or the pressure fluctuation amplitude exceeds the preset pressure threshold, mark the sub-model as a sub-model to be corrected.

[0155] Calculate the deviation rate between the predicted value of the rock mass strength and the designed value of the rock mass strength stored in the sub-model. The designed value of the rock mass strength stored in the sub-model is determined according to geological exploration and engineering requirements during the project design stage. Assume that the designed value of the rock mass strength stored in the sub-model is 100 MPa, and the predicted value of the rock mass strength of the current excavation face is 120 MPa. The calculation method of the deviation rate is to first calculate the difference between the two, that is, 120 MPa - 100 MPa = 20 MPa, and then divide this difference by the designed value of the rock mass strength to obtain a deviation rate of 20 MPa ÷ 100 MPa = 0.2, that is, 20%.

[0156] At the same time, compare the pressure fluctuation amplitude with the preset pressure threshold. The preset pressure threshold is set according to the performance and safety requirements of the tunneling equipment. Assume that the preset pressure threshold is 2.5 MPa, and the current pressure fluctuation amplitude is 3 MPa. Since 3 MPa is greater than 2.5 MPa, it indicates that the pressure fluctuation amplitude exceeds the preset pressure threshold.

[0157] If the deviation rate exceeds the preset deviation rate threshold (assumed to be 15%) or the pressure fluctuation amplitude exceeds the preset pressure threshold, mark the sub-model as a sub-model to be corrected. In this example, the deviation rate of 20% exceeds the preset deviation rate threshold of 15%, and the pressure fluctuation amplitude of 3 MPa exceeds the preset pressure threshold of 2.5 MPa. Therefore, the sub-model is marked as a sub-model to be corrected.

[0158] Step S214: According to the marking instruction of the sub-model to be corrected, arrange dense monitoring points within the spatial range corresponding to the sub-model, and collect data on the joint density of the rock mass and the groundwater seepage rate.

[0159] After the sub-model is marked as a model to be corrected, encrypted monitoring points will be arranged within the spatial range corresponding to the sub-model according to the marking instructions. First, the spatial range of the sub-model is accurately defined to clarify the specific area where encrypted monitoring is required. For example, if the sub-model is a cuboid space with a length of 20 meters, a width of 15 meters, and a height of 10 meters, the positions of the monitoring points are determined based on this cuboid space.

[0160] For the collection of rock joint density data, monitoring points are set within this space in a certain grid layout. Assuming that the points are arranged at a density of one point every 2 meters, then 10 points can be set in the length direction, 8 points in the width direction, and 5 points in the height direction, resulting in a total of 10×8×5 = 400 monitoring points. At each monitoring point, joint measuring tools such as joint compasses are used to carefully measure information such as the number, strike, and dip of the rock joints. By counting the number of joints within a certain range, the rock joint density is calculated. For example, in a cube-shaped rock mass area with a side length of 1 meter, if it is measured that there are 10 joints, then the rock joint density in this area is 10 joints per cubic meter.

[0161] For the collection of groundwater seepage rate data, a piezometer and a flow measurement device are installed at the arranged monitoring points. The piezometer is used to measure the pressure of the groundwater, and the flow measurement device is used to measure the flow rate of the groundwater. By continuously recording the above data over a period of time and combining parameters such as the porosity of the rock mass, the groundwater seepage rate is calculated. For example, at a monitoring point, the piezometer shows that the groundwater pressure is 0.1 MPa, and the flow measurement device records that the amount of water passing through a specific cross-section within 1 hour is 0.5 cubic meters. Based on information such as the porosity of the rock mass and after a series of calculations (considering factors such as the cross-sectional area of the water flow channel and the pressure difference), the groundwater seepage rate at this monitoring point is obtained as 0.1 cubic meters per (hour·square meter).

[0162] Step S215: Integrate the joint density data and seepage rate data collected at the encrypted monitoring points with the fracture development density obtained by inversion to generate updated rock layer gradient change characteristic data.

[0163] Integrate the rock joint density data and groundwater seepage rate data collected at the encrypted monitoring points with the fracture development density data obtained by inversion before. First, preprocess the above data to ensure that they have the same spatial reference and time reference. For example, unify the data collected at different monitoring points at different times into the same time interval and spatial coordinate system.

[0164] For the joint density data and fracture development density data of rock masses, both of which reflect the degree of fragmentation of rock masses, they can be fused by the method of weighted average. Suppose the weight of the joint density data of rock masses is 0.6, and the weight of the fracture development density data is 0.4. In a certain monitoring area, the measured value of the joint density of the rock mass is 8 joints per cubic meter, and the inverted value of the fracture development density is 6 joints per cubic meter. Then the fused fragmentation degree index is 8×0.6 + 6×0.4 = 7.2 joints per cubic meter.

[0165] For the groundwater seepage rate data, which is related to the permeability of rock masses, it can be combined with the fused fragmentation degree index, considering the influence of the fragmentation degree of rock masses on the groundwater seepage rate. For example, by establishing an empirical model, according to the fused fragmentation degree index and the measured groundwater seepage rate, a comprehensive permeability index is calculated. Suppose in a certain area, the fused fragmentation degree index is 7.2 joints per cubic meter, and the groundwater seepage rate is 0.1 cubic meters per (hour·square meter). According to the empirical model, the comprehensive permeability index of this area is calculated to be 0.12 cubic meters per (hour·square meter) (this empirical model may consider factors such as the non-linear relationship between the fragmentation degree and permeability).

[0166] By performing such processing on each monitoring area and integrating all the data, updated data on the gradient change characteristics of rock strata are generated, which includes the spatial distribution of the fused fragmentation degree index, the comprehensive permeability index, etc., reflecting the gradient change characteristics of rock strata at different positions.

[0167] Step S216: Based on the updated data on the gradient change characteristics of rock strata, re-execute the process of dynamic identification of surrounding rock types and generation of safety warning signals. If the updated comprehensive risk value is lower than the original comprehensive risk value, the model correction is completed; otherwise, the step of arranging encrypted monitoring points is re-triggered.

[0168] Based on the updated data on the gradient change characteristics of rock strata, re-execute the dynamic identification process of surrounding rock types. First, input the updated data on the gradient change characteristics of rock strata into a pre-trained surrounding rock classification model. This surrounding rock classification model has been trained with a large amount of historical data and can identify the types of surrounding rock according to various characteristics of the rock strata. After inputting the data, the surrounding rock classification model performs a series of operations, including convolution, pooling, fully connected, etc., and outputs a vector of the probability distribution of surrounding rock categories. For example, the output vector is [0.2, 0.3, 0.5], indicating that the probability of the surrounding rock in this area belonging to the first type is 0.2, the probability of belonging to the second type is 0.3, and the probability of belonging to the third type is 0.5. Select the type with the highest probability as the label of the surrounding rock category in this area.

[0169] Next, according to the new surrounding rock identification results, recalculate the geological risk index of the area where the preset protected object is located. The calculation process is similar to the previous one. Based on the stability coefficient corresponding to the surrounding rock category label, the groundwater pressure value and the in-situ stress concentration degree in the geological parameter set of the sub-model (including updated parameters such as the degree of rock mass fragmentation and permeability), a weighted evaluation model is constructed. Assume that the weight of the stability coefficient is 0.4, the weight of the groundwater pressure value is 0.3, and the weight of the in-situ stress concentration degree is 0.3. During the calculation process, obtain the updated groundwater pressure value and in-situ stress concentration degree data, and combine with the stability coefficient corresponding to the new surrounding rock category label to calculate the updated comprehensive risk value.

[0170] Compare the updated comprehensive risk value with the original comprehensive risk value. If the updated comprehensive risk value is lower than the original comprehensive risk value, it indicates that through densifying monitoring and data updating, the correction of the model has played a positive role, making the assessment of geological risks more accurate. At this time, the model correction is completed. For example, if the original comprehensive risk value is 0.2 and the updated comprehensive risk value is 0.15, it shows that the model has been optimized.

[0171] If the updated comprehensive risk value is not lower than the original comprehensive risk value, it means that the current data updating and model correction are not perfect enough, and there may be some geological conditions that have not been accurately monitored. At this time, re-trigger the step of arranging densified monitoring points. On the basis of the original densified monitoring points, further increase the density of monitoring points or expand the monitoring range. For example, increase the density of monitoring points from one every 2 meters to one every 1 meter, or expand the monitoring range from a partial area of the sub-model to the entire sub-model area, and collect data such as rock mass joint density data and groundwater seepage rate data again, and repeat the subsequent data fusion, surrounding rock type identification, and risk assessment steps until the updated comprehensive risk value is lower than the original comprehensive risk value, and the model correction is completed.

[0172] Figure 2 FIG. shows a schematic diagram of exemplary hardware and software components of a potential hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns provided by some embodiments of the present invention that can implement the idea of the present invention. For example, the processor 120 can be used on the potential hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns and is used to execute the functions in the present invention.

[0173] The potential hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns can be a general-purpose server or a special-purpose server, both of which can be used to implement the method for analyzing potential hazards based on 3D geological intelligent modeling of underground caverns of the present invention. Although only one server is shown in the present invention, for convenience, the functions described in the present invention can be implemented in a distributed manner on multiple similar platforms to balance the processing load.

[0174] For example, a safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns may include a network port 110 connected to a network, one or more processors 120 for executing program instructions, a communication bus 130, and different forms of storage media 140, such as disks, ROM, or RAM, or any combination thereof. Exemplarily, the safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns may also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The method of the present invention can be implemented according to the above program instructions. The safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns further includes an input / output (I / O) interface 150 between the computer and other input / output devices.

[0175] For ease of explanation, only one processor is described in the safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns. However, it should be noted that the safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns in the present invention may also include multiple processors. Therefore, the steps performed by one processor described in the present invention can also be jointly executed or separately executed by multiple processors. For example, if the processor of the safety hazard analysis system 100 based on 3D geological intelligent modeling of underground caverns executes step A and step B, it should be understood that step A and step B can also be jointly executed by two different processors or separately executed in one processor. For example, the first processor executes step A, the second processor executes step B, or the first processor and the second processor jointly execute steps A and B.

[0176] In addition, an embodiment of the present invention further provides a readable storage medium, in which computer-executable instructions are preset. When the processor executes the computer-executable instructions, the safety hazard analysis method based on 3D geological intelligent modeling of underground caverns as described above is implemented.

[0177] It should be noted that, in order to simplify the presentation of the disclosure of the present invention and thus help the understanding of one or more embodiments of the invention, in the foregoing description of the embodiments of the present invention, sometimes multiple features are merged into one embodiment, drawing, or description thereof.

Claims

1. A method for analyzing potential safety hazards based on 3D geological intelligent modeling of underground caverns, characterized in that, The method includes: Based on the three-dimensional model data of the underground chamber, the three-dimensional model is cut into multiple sub-models according to the geological condition segmentation rules, and the segmentation rules dynamically adjust the cutting spacing of each sub-model according to the characteristics of the rock formation gradient change; Perform gradient analysis on the longitudinal rock formation attributes and axial rock formation attributes of each sub-model to generate the rock formation transition characteristic data of the sub-model, and the rock formation transition characteristic data includes the attribute change rate between adjacent rock formations and the spatial distribution of the transition region; According to the rock formation transition characteristic data, dynamically identify the surrounding rock type of each sub-model, and map and match the surrounding rock identification result with the spatial coordinates of the preset protected object to obtain a mapping and matching relationship, and collect the geological parameter set of the sub-model in real time based on the mapping and matching relationship; Trigger a safety warning signal according to the numerical interval of the geological parameter set, and generate an emergency response plan set associated with the safety warning signal.

2. The safety hazard analysis method based on 3D geological intelligent modeling of underground caverns according to claim 1, characterized in that The step of cutting the three-dimensional model into multiple sub-models according to the geological condition segmentation rules based on the three-dimensional model data of the underground chamber includes: Fuse the geological radar detection data, the core drilling data and the UAV oblique photography data to generate an initial three-dimensional geological model. Among them, the geological radar detection data analyzes the distribution of rock mass fissures through the reflection signals of multi-band electromagnetic waves, the core drilling data calibrates the rock formation attributes through the physical and mechanical parameters of the core samples, and the UAV oblique photography data reconstructs the surface and shallow geological structures through multi-view images; Perform spatial grid division on the initial three-dimensional geological model, and screen out the non-uniform regions where the difference in rock mass strength and the difference in permeability coefficient between adjacent grid units exceed the preset threshold based on the rock mass strength distribution curve and the permeability coefficient change rate of the grid unit; Perform secondary data collection on the non-uniform region to obtain supplementary exploration data on the occurrence of rock mass joints and the occurrence state of groundwater, and perform point cloud modeling on the excavation surface through a three-dimensional laser scanner, and spatially superimpose the supplementary exploration data and the point cloud modeling data on the initial three-dimensional geological model to generate corrected three-dimensional model data; According to the characteristics of the rock formation gradient change in the corrected three-dimensional model data, cut the sub-models at a dynamic spacing along the axis direction of the chamber. The adjustment logic of the cutting spacing is: when it is detected that the change rate of the rock mass strength in the adjacent region exceeds the set gradient threshold, reduce the cutting spacing to improve the geological attribute analysis accuracy of the sub-model, and the included angle between the boundary of each sub-model and the geological structure line is kept within a set range through spatial vector calculation to avoid the cutting surface forming an acute intersection with the rock formation strike.

3. The safety hazard analysis method based on 3D geological intelligent modeling of underground caverns according to claim 1, wherein The step of performing gradient analysis on the longitudinal rock formation attributes and axial rock formation attributes of each sub-model to generate the rock formation transition characteristic data of the sub-model includes: Extract the rock mass quality index sequence along the axis direction of the chamber. The rock mass quality index sequence is generated by comprehensively calculating the borehole wave velocity test, the rock mass integrity coefficient and the RQD value. Perform first-order difference calculation on the rock mass quality index sequence to obtain the quality index change gradient between adjacent rock formations, and mark the section where the absolute value of the gradient of the quality index change gradient exceeds the set threshold as the mutation region; Extract the distribution characteristics of the elastic modulus of the rock mass in the direction perpendicular to the axis of the cavern. Through the finite element interpolation algorithm, convert the discrete borehole elastic modulus data in the distribution characteristics of the rock mass elastic modulus into a continuous spatial distribution matrix, calculate the Gaussian curvature of the spatial distribution matrix, and generate a modulus change rate distribution map reflecting the spatial change rate of the rock mass stiffness; Perform a spatial overlay analysis on the position coordinates of the mutation region and the modulus change rate distribution map to determine the main gradient direction of the attribute change rate and the geometric shape of the transition region in the strata transition characteristic data. Among them, the geometric shape is extracted as a polygon patch through the boundary tracking algorithm, and its area, perimeter, and relative position relationship with the boundary of the sub-model are calculated.

4. The safety hazard analysis method based on 3D geological intelligent modeling of underground chambers according to claim 1, characterized in that, Dynamically identify the surrounding rock types of each sub-model according to the strata transition characteristic data, including: Input the strata transition characteristic data into a pre-trained surrounding rock classification model, and output a probability distribution vector of surrounding rock categories. Among them, the surrounding rock classification model is constructed based on the convolutional neural network architecture, and its training data includes sample strata transition characteristic data of different surrounding rock types in historical projects; Generate a surrounding rock identification result with surrounding rock category labels according to the probability distribution vector of surrounding rock categories, and generate a corresponding three-dimensional surrounding rock type distribution heat map in combination with the spatial position coordinates of the sub-model. The rendering rule of the three-dimensional surrounding rock type distribution heat map is: assign set color codes to different surrounding rock categories, and adjust the transparency according to the probability value of the implemented surrounding rock category label to reflect the classification confidence; Optimize the boundary of the region where the difference between the surrounding rock category labels of adjacent sub-models in the surrounding rock type distribution heat map exceeds the set value. Among them, the method of boundary optimization includes: extracting the strata transition characteristic data of the boundary region for secondary classification. If the secondary classification result is consistent with the surrounding rock identification result, the original boundary is maintained; otherwise, the sub-model cutting line is adjusted based on the secondary classification result to generate a corrected surrounding rock identification result.

5. The safety hazard analysis method based on 3D geological intelligent modeling of underground chambers according to claim 1, characterized in that Map and match the surrounding rock identification result with the spatial coordinates of the preset protected object to obtain a mapping and matching relationship, and collect the set of geological parameters of the sub-model in real time based on the mapping and matching relationship, including: Determine the corresponding three-dimensional influence range according to the geographical coordinates of the preset protected object. The three-dimensional influence range is calculated and generated by the geological disturbance propagation model. The input of the implemented geological disturbance propagation model includes the structural type, foundation burial depth, and rock mass stress distribution parameters of the protected object; Project the three-dimensional influence range onto the three-dimensional model, determine the set of sub-models intersecting with it through spatial topological relationship analysis, and establish a mapping and matching relationship between the preset protected object and the sub-model; Collect the set of geological parameters of the sub-model in real time based on the mapping and matching relationship, and extract the surrounding rock identification results of the set of sub-models, and calculate the geological risk index of the area where the preset protected object is located. Among them, the specific calculation logic is: construct a weighted evaluation model based on the stability coefficient corresponding to the surrounding rock category label, the groundwater pressure value and the in-situ stress concentration degree in the set of geological parameters of the sub-model, and output the risk index level; When the risk index level exceeds a preset threshold, a dynamic monitoring link is established between the preset protected object and the corresponding sub-model. The dynamic monitoring link collects real-time data on changes in geological parameters through the deployment of fiber optic strain sensors, piezometers, and microseismic monitoring equipment, and transmits the data on changes in geological parameters to the central processing system for early warning analysis.

6. The safety hazard analysis method based on three-dimensional geological intelligent modeling of underground chambers according to claim 1, wherein, Trigger a safety warning signal according to the numerical range of the geological parameter set, and generate a set of emergency response plans associated with the safety warning signal, including: Normalize the actual measured value of each geological parameter in the geological parameter set with its corresponding safety threshold to generate a deviation coefficient for each parameter; Perform weighted calculation on the deviation coefficients based on a preset deviation weight matrix to generate a comprehensive risk value; When the comprehensive risk value exceeds the first threshold, trigger a first warning signal and generate a preliminary set of response plans, which includes the grouting reinforcement range, support parameter adjustment instructions, and monitoring frequency increase parameters; When the comprehensive risk value exceeds the second threshold, trigger a second warning signal and activate an emergency response process, which includes a stop construction instruction, a drainage system start instruction, and a multi-objective optimization plan generation instruction; According to the first warning signal or the second warning signal, screen the response cases matching the current geological parameter set from the historical database, and extract the plans meeting the preset safety indicators in the response cases as a set of candidate plans; Perform a construction condition matching verification on the set of candidate plans, and eliminate infeasible plans according to the available status of the current construction equipment, the material inventory, and the construction period constraint conditions; Sort the verified set of candidate plans according to the implementation efficiency, cost consumption, and risk suppression effect to generate a final set of emergency response plans.

7. The safety hazard analysis method based on 3D geological intelligent modeling of underground caverns according to claim 1, characterized in that, The method further includes a dynamic update step for the sub-model: Collect real-time vibration spectrum data and hydraulic propulsion pressure data of the tunneling equipment, and extract the main vibration frequency characteristics and the pressure fluctuation amplitude; Input the main vibration frequency characteristics and the pressure fluctuation amplitude into a preset rock mass strength - geological parameter mapping model to invert the predicted value of the rock mass strength and the fracture development density of the current excavation face; Calculate the deviation rate between the predicted value of the rock mass strength and the designed value of the rock mass strength stored in the sub-model, and at the same time compare the pressure fluctuation amplitude with a preset pressure threshold. If the deviation rate or the pressure fluctuation amplitude exceeds the preset pressure threshold, mark the sub-model as a sub-model to be corrected; According to the marking instruction of the sub-model to be corrected, arrange dense monitoring points within the spatial range corresponding to the sub-model, and collect data on the joint density of the rock mass and the groundwater seepage rate; Fuse the joint density data and seepage rate data collected by the dense monitoring points with the inverted fracture development density to generate updated data on the gradient change characteristics of the rock stratum; Based on the updated data on the gradient change characteristics of the rock stratum, re-execute the process of dynamic identification of the surrounding rock type and generation of safety warning signals. If the updated comprehensive risk value is lower than the original comprehensive risk value, the model correction is completed; otherwise, the step of arranging dense monitoring points is re-triggered.

8. The safety hazard analysis method based on 3D geological intelligent modeling of underground chambers according to claim 7, characterized in that, Inputting the vibration main frequency characteristics and the pressure fluctuation amplitude into a preset mapping model of rock mass strength - geological parameters to invert the predicted value of the rock mass strength and the fracture development density of the current excavation face includes: Performing normalization processing on the vibration main frequency characteristics and the pressure fluctuation amplitude respectively to generate a normalized vibration main frequency characteristic and a normalized pressure fluctuation amplitude, wherein the scaling coefficient of the normalization processing is dynamically adjusted according to the range of historical vibration spectrum data and the hydraulic propulsion pressure; Aligning and splicing the normalized vibration main frequency characteristic and the normalized pressure fluctuation amplitude according to the time stamp into a multi - dimensional feature vector, inputting the multi - dimensional feature vector into the rock mass strength - geological parameter mapping model for forward calculation, and outputting the predicted value of the rock mass strength and the fracture development density of the current excavation face, wherein the forward calculation includes mapping the multi - dimensional feature vector to the regression values of the rock mass strength and the fracture density through a fully - connected neural network layer.

9. The safety hazard analysis method based on 3D geological intelligent modeling of underground caverns according to claim 8, characterized in that, The mapping of the multi - dimensional feature vector to the regression values of the rock mass strength and the fracture density through the fully - connected neural network layer includes: Extracting the distribution weight of the vibration main frequency band in the multi - dimensional feature vector and the time - series change slope of the pressure fluctuation amplitude to generate a frequency - domain - time - domain correlation feature matrix; Inputting the frequency - domain - time - domain correlation feature matrix into the hidden units of the fully - connected neural network layer and calculating the output vector of the hidden layer through an activation function; Inputting the output vector of the hidden layer into the rock mass strength regression branch and the fracture density regression branch respectively, wherein the weight parameters of the rock mass strength regression branch are supervised and trained through historical borehole core strength data, and the weight parameters of the fracture density regression branch are supervised and trained through historical borehole fracture imaging data, and outputting the predicted value of the rock mass strength and the fracture development density of the current excavation face.

10. A safety hazard analysis system based on three-dimensional geological intelligent modeling of underground chambers, characterized in that, Including a processor and a memory, the memory is connected to the processor, the memory is used to store programs, instructions or codes, and the processor is used to execute the programs, instructions or codes in the memory to implement the safety hazard analysis method based on three - dimensional geological intelligent modeling of underground caverns according to any one of claims 1 - 9 above.

Images