A method and system for inspecting underground caverns

Abstract

The present application relates to the technical field of route planning, and particularly relates to a method and system for underground cavern inspection. The method comprises the following steps: obtaining geological survey data, microseismic monitoring data and historical deformation monitoring data; performing rock mass viscoelasticity constitutive modeling according to the geological survey data and the microseismic monitoring data to obtain a rock mass viscoelasticity model; performing deformation trend prediction according to the rock mass viscoelasticity model and the historical deformation monitoring data to obtain deformation trend data; performing deformation disturbance resonance analysis according to the deformation trend data and the microseismic monitoring data to obtain deformation disturbance resonance data; performing risk region grading on the deformation disturbance resonance data to obtain risk region grading data; and generating an inspection map according to the risk region grading data to obtain inspection map data, so as to assist in underground cavern inspection operations. The present application can accurately identify potential risk regions, dynamically plan inspection tasks and improve inspection efficiency and decision timeliness.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of route planning, in particular to an underground cavern inspection method and system. BACKGROUND

[0002] With the deepening of the development and utilization of underground space, the long-term stability and safety of structures such as tunnels, mines, and water caverns are facing significant challenges. Especially in complex geological environments, rock mass structures are easily affected by multiple factors such as ground stress disturbance, surrounding rock creep, and joint softening, which can cause abnormal deformation and even catastrophic instability. Therefore, how to accurately identify potential risk areas, dynamically plan inspection tasks, and improve inspection efficiency and decision-making timeliness has become a key problem in the intelligent operation and maintenance of underground engineering. SUMMARY

[0003] To solve the above technical problems, the present application provides an underground cavern inspection method and system to solve at least one of the above technical problems.

[0004] The present application provides an underground cavern inspection method, which comprises:

[0005] S1, obtaining geological survey data, microseismic monitoring data, and historical deformation monitoring data; building a rock mass viscoelastic constitutive model according to the geological survey data and the microseismic monitoring data, and obtaining the rock mass viscoelastic model;

[0006] S2, predicting the deformation trend according to the rock mass viscoelastic model and the historical deformation monitoring data, and obtaining deformation trend data; performing deformation disturbance resonance analysis according to the deformation trend data and the microseismic monitoring data, and obtaining deformation disturbance resonance data;

[0007] S3, classifying the risk areas according to the deformation disturbance resonance data, and obtaining risk area classification data;

[0008] S4, generating an inspection map according to the risk area classification data, and obtaining inspection map data to assist in underground cavern inspection operations.

[0009] In the present application, by fusing geological survey data, microseismic monitoring data, and historical deformation monitoring data, a viscoelastic constitutive model is constructed to realize physical modeling of the response of rock mass structure. Based on the model output and historical data, deformation trend prediction is carried out, and disturbance resonance analysis is carried out in combination with microseismic characteristics, which can identify potential resonance areas and local instability risks in advance. Through risk area classification and intelligent inspection map generation driven by resonance results, the inspection path is more focused on high-risk areas, improving the pertinence and efficiency of the inspection, and having good predictability, real-time performance, and spatial coupling perception ability, which can effectively assist in underground cavern operation and management and disaster warning response.

[0010] Optionally, the S1 comprises:

[0011] obtain geological survey data, microseismic monitoring data and historical deformation monitoring data;

[0012] According to the geological survey data and the microseismic monitoring data, the constitutive model is constructed to obtain the constitutive model;

[0013] According to the microseismic monitoring data, the microseismic feature extraction is carried out to obtain the microseismic feature data;

[0014] The microseismic feature data is processed by the viscous field inversion figure to obtain the viscous field inversion figure data;

[0015] According to the viscous field inversion figure data, the constitutive model is synchronously updated to obtain the fracture network evolution model;

[0016] The fracture network evolution model and the geological survey data are coupled to obtain the rock mass viscoelastic model.

[0017] In the present application, the dynamic disturbance perception mechanism and the fracture evolution feedback path are introduced on the basis of the traditional constitutive model construction, which can more truly simulate the nonlinear mechanical response process of the rock mass after being disturbed in the complex underground environment. Through feature extraction of microseismic monitoring data, including source location distribution, frequency energy ratio and energy release density, etc. The stress concentration and microcrack initiation area in the rock mass can be described in real time. Combined with the regional viscous response atlas constructed by viscous field inversion, the local intensity weakening trend of the microseismic active area is reflected. The constitutive model parameters are locally corrected and time-sequentially updated to generate a fracture network evolution model with time evolution and spatial heterogeneity.

[0018] Optionally, the deformation trend prediction comprises:

[0019] According to the rock mass viscoelastic model and the historical deformation monitoring data, the deformation network structure is constructed to obtain the deformation network representation data;

[0020] The trend change boundary of the deformation network representation data is identified to obtain the trend boundary data;

[0021] The variable response degree of the trend boundary data is evaluated to obtain the trend response degree data;

[0022] According to the trend response degree data, the rock mass viscoelastic model is mapped to obtain the deformation trend number.

[0023] In the present application, by constructing a deformation representation network structure based on the viscoelastic model of rock mass and the fusion of historical monitoring data, the strain evolution characteristics of rock mass under multi-scale and multi-stage are effectively captured. Through trend change boundary identification, the displacement mutation or strain acceleration area can be accurately located, and the identification ability of local unstable evolution is improved. The system performs variable response degree evaluation, quantifies the influence intensity of each disturbance factor on the trend change, and helps to identify the main control factor and abnormal response path. Through the trend mapping mechanism, the analysis results are written back to the viscoelastic model, realizing the spatial prediction and time evolution simulation of future deformation distribution.

[0024] Optionally, the deformation disturbance resonance analysis comprises:

[0025] According to the deformation trend data and the microseismic monitoring data, disturbance same frequency feature extraction is performed to obtain disturbance same frequency feature data;

[0026] According to the disturbance same frequency feature data, disturbance resonance identification is performed to obtain disturbance resonance data;

[0027] Resonance candidate unit extraction is performed on the disturbance resonance data to obtain resonance candidate unit data;

[0028] The resonance candidate unit data is subjected to clustering processing to obtain resonance body data;

[0029] The resonance body data is subjected to resonance intensity grading to obtain resonance grading data;

[0030] According to the resonance grading data, the deformation trend data is mapped to obtain deformation disturbance resonance data.

[0031] In the present application, by extracting disturbance same frequency features from deformation trend data and microseismic monitoring data, the frequency domain coupling relationship between structural response and microseismic disturbance can be effectively identified, overcoming the limitations of traditional methods that rely only on a single physical quantity to determine risks. Through resonance identification and candidate unit extraction, potential resonance areas can be accurately located in space, and clustering algorithms are used to form resonance bodies with structural continuity. Further, by resonance intensity grading, a disturbance level system is constructed, which helps to quantitatively evaluate the instability risk caused by local resonance. By mapping the resonance level results to the deformation trend field, the physical enhancement and spatial feature reconstruction of the trend results are realized, improving the identification ability and forward-looking judgment level of the whole system for potential disaster areas of rock mass.

[0032] Optionally, the disturbance same frequency feature extraction comprises:

[0033] According to the deformation trend data and the microseismic monitoring data, empirical mode decomposition is performed to obtain mode decomposition data;

[0034] Frequency feature extraction is performed on the mode decomposition data to obtain frequency feature data;

[0035] The frequency characteristic data is time-synchronized to obtain frequency characteristic synchronization data;

[0036] The frequency characteristic synchronization data is subjected to main frequency overlap calculation and cross-spectrum coherence extraction to obtain main frequency overlap data and cross-spectrum coherence data, respectively;

[0037] The deformation trend data and the microseismic monitoring data are subjected to coupling degree map construction according to the main frequency overlap data and the cross-spectrum coherence data to obtain coupling degree map data;

[0038] The coupling degree map data is subjected to same-frequency response extraction to obtain disturbance same-frequency characteristic data.

[0039] In the present application, the multi-scale modal characteristics in the nonlinear and non-stationary signals are effectively captured by performing empirical mode decomposition on the deformation trend data and the microseismic monitoring data, and the frequency domain analysis capability for local complex deformation and microseismic events is improved. Through frequency characteristic extraction and time synchronization processing, the different source heterogeneous data are uniformly aligned in the frequency domain dimension, and the comparability of the same-frequency dynamic changes is enhanced. The main frequency overlap and the cross-spectrum coherence indexes are used to accurately depict the synchronism and coupling strength of the deformation and microseismic signals in the key frequency band. The coupling degree map is constructed based on the above indexes, and the high-coupling region is extracted, so that the disturbance same-frequency response can be cooperatively identified in the time-frequency-space dimension.

[0040] Optionally, the disturbance resonance identification comprises:

[0041] The microseismic energy release density data and the deformation gradient intensity data are obtained by performing microseismic energy release density calculation and deformation gradient intensity calculation according to the disturbance same-frequency characteristic data;

[0042] The resonance point data are obtained by performing lightweight resonance judgment according to the microseismic energy release density data and the deformation gradient intensity data;

[0043] The disturbance resonance data are obtained by performing disturbance resonance data generation according to the resonance point data.

[0044] In the present application, the microseismic energy release density and the deformation gradient intensity are extracted from the disturbance same-frequency characteristic data, and the dual representation of microcosmic dynamic disturbance and macroscopic structure response is realized. The microseismic energy release density is used to measure the disturbance intensity of the local rock mass, and the deformation gradient intensity reflects the strain concentration trend in the rock mass, which together constitute an important physical basis for resonance judgment. Through the lightweight resonance judgment mechanism, the high-probability resonance point can be identified under the premise of ensuring real-time, and the resonance identification efficiency and response speed are improved. The generated disturbance resonance data have clear physical meaning and spatial positioning characteristics, and effectively enhance the forward-looking perception ability of the system to the structure resonance instability.

[0045] Optionally, the resonance candidate unit extraction includes:

[0046] The spatial continuity grid is constructed for the disturbance resonance data to obtain resonance point network data;

[0047] The local response screening is performed on the resonance point network data to obtain response hotspot kernel data;

[0048] The spatial topology filtering is performed according to the response hotspot kernel data to obtain resonance candidate unit data.

[0049] In the present application, the spatial continuity grid is constructed for the disturbance resonance data, which can effectively convert the discrete resonance point data into the resonance point network with spatial topology structure, and improve the organization and computability of the data in spatial analysis. The high resonance activity area is extracted through the local response screening, which can quickly focus on the abnormal response hotspot with engineering significance and avoid excessive response to the background disturbance. The system performs spatial topology filtering, judges the connectivity and structural consistency of the response hotspot, removes the isolated noise points and non-structural abnormal areas, and generates the resonance candidate unit with continuous shape and stable structure.

[0050] Optionally, S3 includes:

[0051] The risk feature extraction is performed on the deformation disturbance resonance data to obtain risk feature data;

[0052] The regional risk classification is performed according to the risk feature data to obtain regional risk classification data;

[0053] The grade connectivity correction is performed according to the regional risk classification data to obtain risk connectivity graph data;

[0054] The stability layering is performed according to the risk connectivity graph data to obtain risk region grading data.

[0055] In the present application, the risk feature extraction is performed on the deformation disturbance resonance data, and the multi-dimensional physical indexes such as resonance strength, deformation trend and microseismic activity are calculated, which can comprehensively construct the quantitative feature system of regional risk. Through the regional risk classification, the preliminary determination of different risk level regions is realized, and the accuracy and efficiency of high risk area identification are effectively improved. The system corrects the grade connectivity, combines the spatial adjacency and the topological consistency of the risk level, repairs the isolated high risk points or the boundary fuzzy area, and makes the risk expression more stable and coherent. The system performs the stability layering operation, divides the evolution area, the stable area and the transition area according to the risk level fluctuation trend and evolution potential, realizes the expansion from static risk assessment to dynamic evolution perception. The present application can provide the regional grading atlas with clear structure and reasonable logic for inspection scheduling, support optimization and early warning response, and significantly improve the fine and intelligent level of underground cavern monitoring and management.

[0056] Optionally, S4 comprises:

[0057] According to the risk area classification data, a patrol candidate task unit graph is constructed, and patrol candidate graph data is obtained;

[0058] The node sorting is performed on the patrol candidate graph data, and patrol sorting graph data is obtained;

[0059] The three-dimensional path planning construction is performed on the patrol sorting graph data, and three-dimensional path data is obtained;

[0060] The space constraint correction is performed on the three-dimensional path data, and the patrol graph data is obtained, so as to perform the underground cavern patrol auxiliary operation.

[0061] In the present application, the risk area classification data is used to construct the patrol candidate task unit graph, so that the accurate positioning and task unit expression of the high-risk area can be realized, and the basis support for the fine scheduling of the patrol task is provided. Through the node sorting mechanism, the priority of the patrol task is divided in combination with the risk level, the change trend and the historical patrol interval, so that the resource utilization efficiency and the response timeliness are improved. The system combines the three-dimensional path planning construction, fully considers the complexity and accessibility of the underground cavern space structure, and realizes the three-dimensional feasibility design of the patrol route. The system performs the space constraint correction, excludes the impassable area and optimizes the path cost, so that the generated patrol graph has the actual execution and environmental adaptability.

[0062] Optionally, the present application also provides an underground cavern patrol system for performing the underground cavern patrol method as described above, and the underground cavern patrol system comprises:

[0063] The rock mass modeling and parameter inversion module is used for acquiring geological survey data, microseismic monitoring data and historical deformation monitoring data; the rock mass viscoelastic constitutive modeling is performed according to the geological survey data and the microseismic monitoring data, and a rock mass viscoelastic model is obtained;

[0064] The deformation prediction and disturbance resonance analysis module is used for performing the deformation trend prediction according to the rock mass viscoelastic model and the historical deformation monitoring data, and obtaining deformation trend data; the deformation disturbance resonance analysis is performed according to the deformation trend data and the microseismic monitoring data, and deformation disturbance resonance data is obtained;

[0065] The multi-dimensional risk area classification module is used for performing the risk area classification on the deformation disturbance resonance data, and obtaining risk area classification data;

[0066] The intelligent patrol graph generation and task planning module is used for performing the patrol graph generation according to the risk area classification data, and obtaining patrol graph data, so as to perform the underground cavern patrol auxiliary operation.

[0067] The underground cavern inspection method provided by the application fuses geological survey data, microseismic monitoring data and historical deformation monitoring data, constructs a rock mass viscoelastic constitutive model with time-varying characteristics, and provides a physical basis for truly restoring the response behavior of the rock mass in a complex underground environment. In S1, the microseismic feature extraction and viscous field inversion mechanism are introduced to dynamically identify the potential crack evolution area and update the constitutive parameters synchronously, so as to form a crack network evolution model with disturbance feedback capability. In S2, the trend boundary identification and response degree mapping mechanism are combined with the viscoelastic graph structure to carry out deformation trend prediction, realize time sequence modeling of the progressive deformation of the structure, and introduce the frequency domain response information of the microseismic disturbance through the resonance mechanism analysis, thereby enhancing the discrimination ability of the "hysteresis-suddenness" coupled deformation area. In S3, the multi-dimensional risk features of the disturbance resonance area are extracted, including the microseismic energy release density, the local deformation gradient intensity and the trend change boundary stability index, and a preliminary risk classification graph is constructed based on the spatial position and feature similarity. The regional connectivity correction mechanism is used to reconstruct the weight of the spatial continuity and disturbance coupling strength between the high-risk areas, and eliminate the pseudo-partition effect caused by the monitoring blind area or data dispersion. According to the multi-factor superposition scoring system, the stability of each risk area is layered, including the critical unstable area, the sub-stable area and the allowable fluctuation area, so as to realize the hierarchical and fine expression of the risk area. In S4, the three-dimensional inspection task graph is constructed and the path planning is completed combined with the spatial constraint, so that the inspection task has the task priority visibility, the path continuity and the environmental feasibility, and the inspection efficiency, the accuracy and the safety decision-making ability in the underground cavern space are significantly improved. BRIEF DESCRIPTION OF DRAWINGS

[0068] Other features, objects and advantages of the application will become more apparent from the following detailed description of non-limiting embodiments made with reference to the accompanying drawings:

[0069] Figure 1 A step flowchart of an underground cavern inspection method of an embodiment is shown;

[0070] Figure 2 A step flowchart of a rock mass modeling and parameter inversion method of an embodiment is shown;

[0071] Figure 3 A step flowchart of a deformation trend prediction method of an embodiment is shown;

[0072] Figure 4 A step flowchart of a multi-dimensional risk area hierarchical method of an embodiment is shown;

[0073] Figure 5 A step flowchart of an intelligent inspection graph generation and task planning method of an embodiment is shown;

[0074] The objectives, functional characteristics and advantages of the present application will be further described with reference to the embodiments and the accompanying drawings. DETAILED DESCRIPTION

[0075] The technical method of the present application will be described clearly and completely below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative work fall within the scope of the present application.

[0076] In addition, the accompanying drawings are only schematic illustrations of the present application, and are not necessarily drawn to scale. The functional entities can be implemented in software, or in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0077] It should be understood that although the terms "first", "second" and the like can be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, without departing from the scope of the example embodiments, a first element can be called a second element, and similarly, a second element can be called a first element. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0078] Referring to Figures 1 to 5 The present application provides a method for inspecting underground caverns, the method comprising:

[0079] S1, obtaining geological survey data, microseismic monitoring data and historical deformation monitoring data; performing rock mass viscoelasticity constitutive modeling according to the geological survey data and the microseismic monitoring data to obtain a rock mass viscoelasticity model;

[0080] In an embodiment, the geological survey data includes lithology profile information, joint and fracture structure parameters, initial ground stress distribution and underground water level conditions; the microseismic monitoring data is collected in a time sequence event point set format, including event occurrence time, magnitude, spatial coordinates, released energy value and spectral feature information; the historical deformation monitoring data is recorded in a set of time sequence monitoring point sets, and the data fields include monitoring time, measurement point number, displacement value or strain value, and three-dimensional spatial position coordinates of the corresponding measurement point. The system constructs an initial rock mass viscoelasticity constitutive model. The model models the stress response as a superposition of elastic response and viscous response, and the specific expression is: wherein is the elastic modulus, used to represent the elastic response ability of the rock mass to deformation; The viscosity coefficient is used to characterize the damping effect of the rock mass over time. The system performs viscous field inversion driven by microseismic events. Microseismic events are divided into three-dimensional grid cells according to spatial location (e.g., with a resolution of 1 cubic meter), and the energy release density per unit time for each grid is statistically analyzed to construct a spatial energy field distribution map. The system uses this energy density field as the optimization objective, adjusting the spatial viscosity coefficient accordingly. To minimize the error between the deformation trend simulated by the viscoelastic model and historical monitoring data, error optimization can be achieved using the least squares method or Bayesian inversion method. A fracture network evolution model is constructed. The spatial intersection of high energy density areas and areas with dense microseismic events is identified and marked as "rock softening zones." Additional fault zone structural nodes are introduced into these areas, and the connectivity graph structure within the regions is updated, forming a composite model that simultaneously characterizes fracture propagation and viscoelastic response properties.

[0081] S2. Based on the rock mass viscoelastic model and historical deformation monitoring data, deformation trend prediction is performed to obtain deformation trend data; based on the deformation trend data and microseismic monitoring data, deformation disturbance resonance analysis is performed to obtain deformation disturbance resonance data.

[0082] In one embodiment, the system constructs a deformation trend prediction model based on a rock mass viscoelastic model and historical deformation monitoring data. The rock mass region is divided into a spatial grid and treated as nodes in a graph structure. Each node carries constitutive parameters (such as elastic modulus) for its corresponding location. viscosity coefficient ) and its historical deformation response sequence. For each node in the graph Construct feature vectors The viscoelastic parameters of the node and the historical strain values at multiple time points are included, which are used to characterize the local time-varying response behavior of the rock mass at the location. The system uses a time convolution network (TCN) to model the above historical sequence. With the historical displacement sequence of each monitoring point as input, the continuous strain evolution trend in the future period (such as 30 days) is predicted, and the output is the predicted strain sequence. In this way, the spatial distribution trend data of the rock mass deformation in the region within the future time period, i.e., the deformation trend data, can be obtained. The system uses the deformation trend data and microseismic monitoring data to jointly perform deformation disturbance resonance analysis. The energy time series signal of each microseismic event and the predicted strain curve of the adjacent position are respectively subjected to empirical mode decomposition to extract multiple intrinsic mode functions (IMF) for representing the signal fluctuation characteristics at different frequency components. The Hilbert transform is performed on each IMF modal signal to extract its main frequency information and instantaneous energy value to form a frequency-energy pair. If the main frequencies of the two signals fall within a frequency interval of ±0.5 Hz from each other, it is counted as one main frequency overlap, and the proportion of the overlapping frequency band in all modes is counted. The cross-spectral coherence function is used to evaluate the coupling strength of the two signals in the frequency domain, and the score ranges from 0 to 1, and the larger the value, the stronger the frequency synchronization. If the main frequency overlap between the two signals is more than 60%, and the cross-spectral coherence of the corresponding frequency band is greater than or equal to 0.7, it is determined that there is a disturbance resonance phenomenon at the position. The identified resonance points constitute a set of spatial events, each event including position coordinates, frequency coupling value, frequency pair information, and corresponding local energy density.

[0083] S3, classifying the risk area based on the deformation disturbance resonance data to obtain risk area classification data;

[0084] In an embodiment, a feature vector constituting a risk assessment index is extracted for each spatial unit (such as a 3D grid voxel). The vector consists of four dimensions, including a deformation rate value representing the strain or displacement change rate of the unit per unit time; a microseismic energy density value representing the microseismic energy released per unit volume within a certain time window; a resonance score value representing the frequency overlap intensity and coherence level score obtained by the unit in the disturbance resonance determination; and a strain gradient value representing the strain distribution change rate within the unit or between adjacent units. The feature vector is used to construct a risk identification vector wherein is the deformation rate value, is the microseismic energy density, is the resonance score, is the strain gradient. The system classifies and evaluates the risk based on a preset classification rule. The classification adopts a five-level risk grading system, and the specific grading standards are as follows: Level 1 (extremely high risk): when the resonance score is greater than 0.9, the deformation rate is greater than 0.8, and the strain gradient is greater than 0.6 Level 2 (high risk): when the resonance score is greater than 0.8, the deformation rate is greater than 0.6, and the strain gradient is greater than 0.4 higher than 10 mm / day, and microseismic energy density exceeds a set high-risk threshold value; Level 3 (high risk): each index is slightly lower than Level 2; Level 2 (medium risk), Level 1 (low risk), Level 1 (extremely low risk): the judgment threshold values of each index are set in descending order according to the risk level. After completing the preliminary risk classification, the system performs connectivity correction processing on the spatially distributed risk areas. If the Euclidean distance between two spatial units is less than 2 meters and the difference in risk level between the two is not more than 1 level, they are considered to be the same connected risk area and are merged; if a high-risk unit is in an isolated state and the number of connected units is less than 3 voxels, the point will be marked as a "boundary drift point" or automatically downgraded to exclude isolated abnormal points from interfering with the overall evaluation. To distinguish the stability state of the risk area, the system calculates the risk volatility index. This index is based on the historical change data of the risk score in the past week (for example, the past 7 days) to calculate its standard deviation wherein is the risk volatility index, is the standard deviation function, is the risk identification vector, is the time variable. If the volatility of a certain area is large, indicating that its risk state is unstable, it is classified as a "dynamic unstable area"; otherwise, if the volatility is small but the risk level itself is high, it is marked as an "evolutionary high-stability area", i.e. the risk is high but the state is relatively continuous, which is suitable for intensive monitoring.

[0085] S4, generate an inspection graph based on the risk area classification data to obtain inspection graph data for assisting in underground cave inspection.

[0086] In an embodiment, the system performs inspection candidate unit construction. The system selects all risk levels not lower than medium risk (i.e. and is divided into a "dynamic instability zone", it is taken as a potential inspection target area. The selected units constitute the initial node set of the inspection graph. The system performs priority ranking on the inspection nodes. Each node is scored according to the following three factors: first, the risk level score of the unit; second, the spatial distance from the current main passage path; third, the time interval since the last inspection of the node. The system sets the weight factors of the three indicators as 0.5, 0.3, and 0.2 respectively, and calculates the priority score value of each node. The node with a higher score represents a more prominent risk or a long-term un-covered area, which will be preferentially included in the inspection path. The system performs three-dimensional path planning construction. The system performs path planning operations based on the three-dimensional structure grid of the underground cavern (such as the geological model or the BIM model). In specific implementation, the system uses a shortest path search algorithm based on graph structure (such as Dijkstra algorithm) to calculate the optimal inspection path between the start and end nodes in the constructed spatial topology graph. The path planning needs to meet the following constraints, such as the total length of each inspection path segment should not exceed 300 meters; the path should not pass through construction barriers, support isolation areas and other forbidden areas (judged according to the pre-set attribute labels in the BIM model); the system will preferentially select areas passing through deployed sensing devices (such as sensors, monitors) to enhance the information collection density along the path. Perform path smoothing correction and graph structure output. The system performs geometric smoothing processing on the initially generated inspection path to improve the actual passage efficiency of the robot or personnel. For example, a B-spline curve fitting method can be used to reconstruct a continuous curve from a polyline path, eliminating unnecessary sharp corners. The system forms a structured inspection graph data, which includes an ordered sequence of inspection task points, corresponding three-dimensional path coordinate sets, and task type labels (such as photographing, laser measurement, support facility inspection, etc.) of each task point. The output data can be encoded in a common format, such as GeoJSON format, glTF format or XML format, to adapt to the visualization and scheduling interface in mobile terminals, smart helmets or ground control systems.

[0087] Optionally, the S1 comprises:

[0088] S11, acquiring geological survey data, microseismic monitoring data, and historical deformation monitoring data;

[0089] In an embodiment, the geological survey data includes the following types of data content, including lithology distribution information, describing the rock type, layer thickness parameters and the profile number of each rock layer in the underground; joint and fracture structure information, including the strike of the joint (i.e. the angle with the north direction), the dip angle (i.e. the angle between the joint surface and the horizontal plane) and the joint density (the number of joints per unit volume); ground stress monitoring data, including the spatial coordinates of the ground stress measurement points, the principal stress direction and the principal stress value; hydrogeological information, describing the water-bearing property of each rock layer and its corresponding permeability coefficient; the above data can be stored and analyzed through Excel table, JSON structured text or Shapefile format in geographic information system, facilitating sharing and fusion among multi-source platforms. The microseismic monitoring data is used to realize real-time sensing of underground stress changes and fracture expansion, and the structured content includes the unique identification number of each event; the event occurrence timestamp (the time accuracy can reach seconds); the three-dimensional spatial coordinates of the event hypocenter ; the magnitude of the microseismic event (generally less than 2.0); the energy release value (unit: joule); the dominant frequency (unit: hertz), which is used to judge the event properties and propagation characteristics. The historical deformation monitoring data reflects the deformation trend of the underground structure over time, and the data content includes the sampling time corresponding to each monitoring time; the monitoring point number and its spatial coordinates; and the displacement vector corresponding to the time (including the direction and amplitude of displacement), which is used to deduce the internal deformation law of the rock mass; the data structure supports spatial voxel division based on a three-dimensional point cloud model, or a structured table is constructed in the form of a gridded monitoring point to realize continuous monitoring under the time and space resolution.

[0090] S12, constructing a constitutive model according to the geological survey data and the microseismic monitoring data, to obtain the constitutive model;

[0091] In an embodiment, the initial form of the viscoelastic model is established as follows: , wherein is the constitutive model formula, such as a stress function, is the value assigned according to the rock type table (for example, granite = 50 GPa), is a strain function, is the initial viscosity modulus set according to the construction disturbance area (for example, 1-10 GPa·s), is a time variable; the space is mapped to a three-dimensional voxel grid to form an initial constitutive tensor field.

[0092] S13, extracting microseismic characteristics according to the microseismic monitoring data, to obtain microseismic characteristic data;

[0093] In an embodiment, the system takes three-dimensional space voxel (such as a cubic unit with each side length of 1 meter) as the basic unit, and performs statistics and analysis on the characteristics of microseismic events in different spatial positions and time windows. The extracted microseismic characteristics include the following four types, including microseismic energy density calculation. The system calculates the total energy released by all microseismic events in a unit time in each spatial voxel as the microseismic energy density of the region. The energy values of all microseismic events falling into the voxel and occurring in the current time window are summed. The energy value can be used as an approximate indicator of the degree of local stress release. The dominant frequency statistics are performed on the microseismic dominant frequency (i.e. event dominant frequency) in each voxel. The mean value of the dominant frequency reflects the concentrated frequency of the microseismic fluctuation in the region, and the standard deviation reflects the dispersion degree of the frequency spectrum distribution. If the mean value of the dominant frequency of a certain region is significantly higher or the standard deviation is significantly abnormal, it can be marked as a high-frequency abnormal region. The event density field construction is to record the number of microseismic events in a unit time window in each voxel and divide it by the spatial volume of the voxel to obtain the event density index. The index reflects the frequency of microseismic activity in a unit space and time, and can be used to identify microseismic active areas or strong fracture activity areas. The system spatially fuses the above-mentioned characteristics to construct a microseismic characteristic tensor data. Each spatial voxel will be associated with a characteristic vector containing four components, which are microseismic total energy, mean value of dominant frequency, standard deviation of dominant frequency, and microseismic event density. The formed microseismic characteristic tensor is a four-dimensional structure that changes with time.

[0094] S14, perform viscous field inversion map processing on the microseismic characteristic data to obtain viscous field inversion map data;

[0095] In an embodiment, the system takes microseismic energy density as the target value of the viscous coefficient inversion model, that is, on the basis of the observed microseismic energy in each spatial voxel, the corresponding viscous parameter is inversely calculated so that the energy density obtained by simulation calculation can fit the true observation value as much as possible. The system constructs the following viscous field inversion optimization model. In a three-dimensional space grid, a viscous parameter value to be solved is set for each voxel unit, denoted as ; based on the current viscous coefficient field , the theoretical energy density simulated by the constitutive model is ; the error is calculated between the simulation value and the actually observed microseismic energy density ; the objective function (loss function) is defined as the sum of the error squares between the simulation value and the observed value in all voxels. That is, the overall error square sum is minimized: the difference between the simulated microseismic energy density and the actual observation value is minimized. The system solves the above objective function, including particle swarm optimization (PSO); L-BFGS algorithm in quasi-Newton method; Bayesian optimization method. Through the above optimization process, the system can obtain the optimal viscous parameter value at each spatial voxel position, denoted as , and a viscous field inversion map covering the whole area is constructed. The system constructs a regularization term to suppress the sudden change of the viscous parameter between adjacent voxels, to control the overall spatial smoothness of the viscous field, and to improve the extrapolation ability and physical consistency of the model in the unobserved area. The output of the system is a three-dimensional tensor field of the viscous field inversion map.

[0096] S15, synchronously updating the constitutive model according to the viscous field inversion map data to obtain a fracture network evolution model;

[0097] In an embodiment, after updating the viscous parameter field, a microseismic concentration area is introduced as an evolution driving area, and it is judged whether the following fracture evolution conditions are met. The system takes the following two criteria as the trigger conditions for fracture evolution: the microseismic event density (i.e. the number of events per unit space per unit time) exceeds the set threshold, denoted as if the event density is greater than the event density threshold ; the time rate of change of the microseismic energy density (i.e. the energy release rate) exceeds the energy threshold, denoted as if the energy density change rate is greater than the set threshold ; the area meeting the above two conditions is marked by the system as an active fracture core area, which is considered as the starting position of microcrack propagation or rock mass failure in the current stage. The system takes the active fracture core as the starting point, combines the known joint structure direction and fracture strike information in the geological survey data, and constructs the fracture propagation path. The fracture propagation is specifically simulating the propagation of the fracture along the direction of the minimum resistance path; preferentially performing linear or sheet fracture structure deduction in the direction of the joint surface; simulating the expansion starting point, endpoint, path and width change of the fracture when it passes through the voxel grid. The processing result can be output as a list of fracture segments, including the spatial starting and ending points, path segment information and estimated fracture width of each fracture, as part of the edge set of the initial fracture network graph. Based on the above static fracture segment, the system establishes a time-series fracture growth model to support dynamic updating of the expansion, connection and healing behavior of the fracture. The model can deduce the spatial growth trend of the fracture segment in time sequence, and adjust in real time in combination with the viscoelastic response and microseismic event driving. In the viscoelastic graph structure, the formed fracture segment is regarded as a "brittle connection edge" with low viscosity and low bearing capacity, and is coupled with the mechanical response of the adjacent area for modeling, and the system outputs a fracture network evolution graph model.

[0098] S16, coupling graph modeling of the fracture network evolution model and the geological survey data to obtain a rock mass viscoelastic model.

[0099] In an embodiment, the system performs grid processing on the entire underground rock mass region according to a preset three-dimensional space division rule to construct a basic three-dimensional graph structure. The nodes of the graph are constructed as each three-dimensional grid voxel unit is regarded as a node in the graph to represent the viscoelastic response unit at the space position. Each node is attached with the following attribute parameters, such as elastic modulus : derived from geological survey data or previous modeling; viscosity modulus : derived from microseismic driving inversion calculation; fracture core label: indicating whether the voxel is an active fracture core region; lithology label: representing the rock type of the rock layer where the voxel is located, serving as semantic auxiliary information for attribute propagation in the graph structure. The edges of the graph are constructed as the connection relationship between nodes, represented as edges in the graph, and each edge represents a certain type of structural connection or physical coupling path. According to different connection mechanisms, the edges can be divided into fracture connection edges, representing actual fracture paths derived from the fracture network evolution model; joint connection edges, representing potential structural weak plane connections between adjacent nodes according to the joint distribution direction in geological survey; and deformation coupling edges, used to represent the mechanical coupling relationship between adjacent viscoelastic units, based on shared boundaries or adjacent field responses. All edges have directionality (such as fracture propagation direction, joint strike, etc.) and weight (which can be assigned according to fracture width, stress gradient, or structural brittleness index) to quantify the coupling strength and propagation bias between structures.

[0100] Optionally, the deformation trend prediction comprises:

[0101] S21, constructing a deformation characterization network structure according to the rock mass viscoelastic model and historical deformation monitoring data to obtain deformation network characterization data;

[0102] In an embodiment, the system performs graph structure modeling. The three-dimensional rock mass region is divided into regular voxel units or hexahedral grid structures, and each voxel unit is regarded as a node in the graph structure to form a node set. Each node is attached with multiple attribute fields, including but not limited to local viscoelastic constitutive parameters (such as elastic modulus and viscosity coefficient), lithology classification label, and deformation response history data sequence of the corresponding region in the near future. The history data sequence can include displacement values or strain values at multiple times within a time window, constituting a time series feature vector. The connection relationship between nodes in the graph structure is constructed. The standard for establishing edge relationship between nodes is based on the spatial proximity and physical attribute similarity between nodes. If the spatial units represented by two nodes are adjacent, or their physical attributes (such as viscosity parameters) differ within a threshold range, an edge relationship can be established between them. The weight of the edge can be calculated by a combination function of spatial distance and attribute difference, for example, in the form of weighted inverse of distance and viscosity coefficient difference, i.e. , wherein is the weight of the edge, is the node spatial position vector of the node, spatial position vector of the node spatial position vector of the node, spatial position vector of the node viscosity modulus value of the node, viscosity modulus value of the node viscosity modulus value of the node. The known joint surface or fracture structure in the rock mass is obtained as a topological enhancement means, which is regarded as an edge with semantic information, so as to improve the geological reality of the structure representation. Time sequence embedding construction is performed. For the historical deformation response data of each node, the deformation features at multiple moments are extracted based on a preset fixed-length time sequence, and a multi-dimensional time feature tensor is constructed. Optionally, the system performs time embedding, for example, based on a position coding method, to encode the time features into auxiliary vectors. The system outputs a deformation representation network structure, including a constructed graph structure ontology (a set of nodes and edges), multi-moment feature information attached to each node, and potential semantic structure information, forming a deformation representation graph data with spatio-temporal expression capability.

[0103] S22, trend change boundary identification is performed on the deformation network representation data to obtain trend boundary data;

[0104] In an embodiment, the system performs deformation trend curve extraction. The system processes the historical displacement or strain sequence carried by each node in the graph structure using a smoothing filter method, for example, selects a Savitzky-Golay filter to remove high-frequency noise and retain the main trend shape. The filtered curve represents the deformation trend evolution process of the spatial unit within the selected time window. Based on the above smoothed curve, the system calculates the first and second derivatives, which correspond to the deformation rate (i.e., displacement change rate) and deformation acceleration (i.e., trend change rate) of the node, respectively. The first derivative is used to capture the overall change direction and rate, and the second derivative is used to depict the acceleration or turning amplitude of the trend. The system identifies feature points of the deformation trend curve of each node according to the set boundary determination rule. If a time point meets any of the following conditions, it can be determined as a trend change boundary. The deformation acceleration value at the current moment exceeds the preset threshold, indicating that the trend has mutated or accelerated sharply. The product of the current deformation rate and the deformation rate at the previous moment is negative, indicating that the trend direction has reversed, i.e., from growth to contraction or from contraction to growth. Within the specified time window, the deviation of the local change rate from the average change rate exceeds the set proportion (such as plus or minus 30%), indicating that a mutation or local instability has occurred. For the time points that meet the above conditions, the system marks them as trend change boundary points and outputs a set of trend boundary points. Each boundary point can be attached to multiple attribute information, such as "change amplitude", "change type", "occurrence time", etc.

[0105] S23, variable response degree evaluation is performed on the trend boundary data to obtain trend response degree data;

[0106] In an embodiment, the system constructs a set of input features for each spatial node. The selected features cover both mechanical parameters and structural evolution attributes of the rock mass, including but not limited to elastic modulus, viscous coefficient, microseismic event density (representing the concentration of microseismic events in a unit space), fracture indicator (used to indicate whether the node is in the fracture distribution area), and deformation trend value of adjacent upstream nodes in the graph structure (used to capture spatial dependence). The above variables form a high-dimensional feature vector for each node. The system sets a response variable to measure the change amplitude before and after the trend jump. Taking the time point when the trend boundary is located as the center, the difference in deformation rate of the corresponding node is calculated at the adjacent two time points before and after it, as the response strength indicator. This indicator represents the severity of the trend jump. The system uses an interpretive modeling tool to model and analyze the causal relationship between the input features and the response variable, including but not limited to a linear regression model for obtaining the linear contribution coefficient of each variable; a SHAP value decomposition algorithm based on game theory for quantifying the marginal contribution of each variable in model prediction; and a gradient boosting decision tree (GBDT) model for measuring the influence degree of variables through feature importance score. Through the above methods, the system can obtain a set of trend response degree score vectors corresponding to each spatial unit, which is composed of the contribution weights of each variable, for example , respectively representing the relative response degrees of elastic modulus, viscous coefficient, microseismic density, fracture distribution, and adjacent influence. The system outputs the response degree vector of each node and marks the "main control factor" field according to the score results, i.e., the parameter variable that plays a dominant role in the trend jump of the node.

[0107] S24, mapping the trend of the rock mass viscoelastic model according to the trend response degree data to obtain a deformation trend number.

[0108] In an embodiment, the trend prediction result is mapped back to the rock mass structure with the response sensitive features, realizing the enhancement of physical correlation of deformation trend and the expression of spatial field. For the predicted sequence, a response degree weighted correction term is introduced: , where is the predicted deformation value of the th node at time , is the initial predicted deformation trend value of the th node at time , is the response influence factor index (such as the risk source or disturbance factor number), is the weighted coefficient of the th factor (response degree weight), is the response degree of the a response value of a factor (such as hot spot activity, risk level, resonance intensity, etc.); the system maps the prediction result back to the three-dimensional structure of the rock mass, and constructs a continuous trend field through three-dimensional interpolation (such as IDW, Kriging); the system calculates and synchronously outputs a trend principal axis direction diagram (based on gradient analysis), a trend high-value area thermal diagram (based on a preset threshold to highlight areas with higher deformation amplitude or deformation rate), and a trend variable voxel diagram (in each voxel position , a fixed time sliding window is selected , the variance of the deformation trend value in the time period is calculated If the variance exceeds a set threshold , it is considered that the region has trend instability. The first-order difference of the trend time series is calculated , wherein is the first-order difference value of the trend at time , representing the deformation increment between the current time and the previous time, is the original deformation trend value at time , is the original deformation trend value at time , and whether a sudden peak value (for example, a fluctuation more than 2 times the mean value occurs multiple times) appears in the sliding window is analyzed, which can be marked as a mutation area. If a voxel satisfies any one of the conditions (high variance or multiple sudden jumps), it is marked as a high-risk area in the trend variable voxel diagram, to prompt the risk area of future significant structural deformation. The spatial trend number field can be output as NetCDF / VTK / GeoTIFF.

[0109] Optionally, the deformation disturbance resonance analysis comprises:

[0110] According to the deformation trend data and the microseismic monitoring data, disturbance same-frequency feature extraction is performed to obtain disturbance same-frequency feature data;

[0111] In an embodiment, the system respectively performs empirical mode decomposition on the deformation trend data and the microseismic signal data to extract respective intrinsic mode function sequences IMF. For each IMF mode, the system performs frequency analysis to extract the following two key frequency domain indicators, including instantaneous frequency, which is used to describe the dominant frequency of the signal at a specific time point; and instantaneous energy spectrum, which represents the local energy intensity of the signal at the frequency. The system performs pairwise matching analysis on the decomposed IMF modes of the deformation trend signal and the microseismic signal to evaluate their coupling degree in frequency, mainly through the following two frequency domain indicators, including the dominant frequency overlap degree indicator, which calculates the ratio of the difference and the average of the dominant frequencies of the deformation signal and the microseismic signal in the corresponding IMF pair. If the ratio is lower than a preset threshold (for example, 0.2), it is determined that the dominant frequencies of the two signals are highly consistent. The specific calculation expression is that the dominant frequency overlap degree is equal to the ratio of the difference and the average of the dominant frequencies of the two signals. When the dominant frequency overlap degree is less than the set threshold (for example, 0.2), it is marked as “dominant frequency consistent”. The cross-spectral coherence degree indicator is calculated, which is used to measure the energy coupling strength of the two signals at a specific frequency. The cross-spectral density and the respective power spectral density are normalized to obtain a coherence value between 0 and 1. The higher the coherence, the stronger the linear coupling relationship between the two signals at the frequency. Specifically, the cross-spectral density between the deformation signal and the microseismic signal is calculated; the power spectral density of each signal is calculated, that is wherein is the cross-spectral coherence function, is the cross-spectral density, is the deformation signal power spectral density, is the microseismic signal power spectral density; and a normalized frequency coherence function is obtained, which is used to quantitatively identify the strong coupling frequency band. The system marks the frequency band that simultaneously satisfies the “dominant frequency overlap” and “high coherence” conditions as the same frequency coupling region.

[0112] According to the disturbance same frequency feature data, disturbance resonance recognition is performed to obtain disturbance resonance data;

[0113] In an embodiment, the system takes the spatial voxels marked as the same frequency coupling region as the input range, and sequentially calculates the following physical parameters in the range, including the microseismic energy release density (E) , which represents the total energy release amount of microseismic per unit volume and per unit time. The system accumulates the energy of all microseismic events falling into the voxel in the current time window, and calculates the local stress release intensity of the region according to the energy. The deformation gradient intensity (G ) is calculated. The system calculates the local gradient field in the spatial voxel based on the three-dimensional displacement field output by the previous trend prediction module. This indicator can be used to identify deformation mutation regions or displacement concentration regions, and is sensitive to the nonlinear response of the structure. The system classifies the above regions according to the following disturbance resonance judgment rules, including the high energy release condition Exceeding the preset energy threshold If the deformation gradient intensity is high, then the region is considered to have high-intensity microseismic activity; the high deformation gradient condition is if the deformation gradient intensity is high. Greater than the threshold The system indicates that a significant local structural response has occurred in the region; the dominant frequency coupling condition is a dominant frequency overlap of less than or equal to 0.2, indicating a high degree of consistency between the microseismic event and the deformation at the dominant frequency; the cross-spectral coherence condition is a cross-spectral coherence of greater than or equal to 0.7, indicating a strong frequency coupling relationship between the two signals. If a spatial voxel simultaneously meets all four of the above composite conditions, the system marks the voxel as a perturbation resonance point, indicating that there may be microseismic-induced resonance or local instability in the region. The system aggregates all identified resonance points according to their spatial location and timestamp, forming a perturbation resonance event set.

[0114] Resonance candidate elements are extracted from the perturbation resonance data to obtain resonance candidate element data;

[0115] In one embodiment, a 3D spatial mesh structure is constructed (e.g., voxel resolution 1m³); the criterion for determining the clustering of resonance points is the number of resonance points within each voxel > 1. (Preset threshold data for the number of resonance points, ranging from 3 to 5) or the average resonance score > (A preset resonance score average threshold, ranging from 0.6 to 0.8), is then marked as a high-response unit. Spatial connectivity determination: The six / eight-way adjacency principle is used to determine whether adjacent voxels are connected; the volume of the connected component is retained to be ≥ The structure (connected component volume threshold data, with values ​​ranging from 5 to 10) outputs resonant candidate units.

[0116] Clustering is performed on the candidate resonance units to obtain the resonance body data;

[0117] In one embodiment, candidate units that are spatially close and have similar features are grouped into a single resonator. Feature vector construction: Location ,in The average microseismic energy density is the average energy of all perturbation resonance points within the candidate element. The average deformation gradient intensity, The dominant frequency mean, The mean cross-spectral coherence is... For the three-dimensional spatial coordinates, use the DBSCAN clustering algorithm (density clustering), with the parameters set to... =2m (radius), MinPts=3; Output connected component clustering results. The obtained resonant data includes spatial center coordinates; average resonant intensity; spatial extension size and volume; dominant frequency standard deviation (representing consistency).

[0118] The resonance intensity of the resonance body data is graded to obtain resonance grading data;

[0119] In an embodiment, a resonance risk level is divided according to a plurality of resonance characteristics. A resonance intensity scoring function is calculated: , wherein is the intensity score value of the resonance body, is an energy term weighting coefficient, and the value is 0.35, is an energy integral index (such as total energy or peak energy) of the resonance body, is a displacement gradient term weighting coefficient, and the value is 0.25, is the displacement gradient modulus value at the boundary of the resonance body, is a connectivity term weighting coefficient, and the value is 0.2, is a topological connectivity index (such as the number of voxels or the size of a subgraph) of the resonance body, is a frequency term weighting coefficient, and the value is 0.2, is the dominant frequency (or main frequency amplitude) of the resonance body; and a level grading rule is set: (very high): > 0.9; (high): 0.75-0.9; (medium): 0.5-0.75; (low): 0.3-0.5; (very low): <0.3; and the grading is marked on the resonance body data, and resonance grading data is obtained.

[0120] According to the resonance grading data, deformation trend data is mapped to obtain deformation disturbance resonance data.

[0121] In an embodiment, the resonance recognition result is reflected to the original deformation trend field. The predicted trend tensor is marked: if the point is located in the resonance body region, an additional field, such as a resonance mark field, is set to “true”, which is used to indicate that the point is in the resonance influence range; a risk level field records the risk level label (such as L5, L4, etc.) of the resonance body corresponding to the point, which reflects the resonance risk intensity of the region. The deformation trend value is adjusted (weighted disturbance): , is a resonance level mapping coefficient, wherein is a deformation trend tensor field after resonance enhancement, is an original deformation trend tensor field, is a disturbance coefficient, which is mapped and set according to the resonance risk level . For example: (very high risk) corresponds to = 0.3; (High risk) Corresponding =0.2; (Medium risk) Corresponding =0.1; (Low risk) Corresponding =0.075; (Very low risk) corresponding =0.05. Output deformation trend data field including resonance enhancement.

[0122] Optionally, the extraction of the perturbation frequency feature includes:

[0123] Empirical mode decomposition was performed based on deformation trend data and microseismic monitoring data to obtain mode decomposition data;

[0124] In one embodiment, the input deformation trend data is a continuous time series. Microseismic monitoring data, aggregated by time period to obtain microseismic event energy sequences. Empirical Mode Decomposition (EMD) is then performed on both time series. For each signal, EMD is performed, decomposing it into several Integral Mode Factors (IMFs). ,in For the deformation trend time series, For the order term of the intrinsic mode functions, The number of intrinsic mode functions (IMFs) obtained from the decomposition. The first deformation trend signal Each intrinsic mode function (IMF) component, This represents the residual term of the deformation trend signal. The modal decomposition of the microseismic signal takes a similar form. Change to .

[0125] Frequency feature data is obtained by extracting frequency features from the mode decomposition data;

[0126] In one embodiment, a Hilbert transform is performed on each IMF to extract the instantaneous frequency. This represents the dominant vibrational frequency and energy of the mode at each moment. The square of the instantaneous amplitude represents the local modal energy intensity; the dominant frequency is calculated as follows: ,in For the first Energy-weighted average frequency of each IMF mode For the first The instantaneous frequency of each IMF mode, For the first The instantaneous energy of each IMF mode, For time variables; construct frequency feature terms for each IMF: Modal energy, duration ,in For the first Energy-weighted average frequency of each IMF mode For the first The frequency standard deviation of each IMF mode For the first The total energy of each IMF mode (which can be...) ), For the first The duration of significant presence of each IMF modal signal.

[0127] Time synchronization is performed on the frequency characteristic data to obtain frequency characteristic synchronization data;

[0128] In one embodiment, the deformation trend and microseismic event signals are resampled to a unified time step (e.g., =1min); Time alignment of the modal dominant frequency sequence: Use a fixed time series (e.g., 60 minutes) to align the modal dominant frequency sequence. and The frequency features are spliced ​​together according to the time series; for cases where the start time of some modes is not synchronized or there are sampling gaps, forward padding or median interpolation is used to align the gaps and obtain synchronized frequency feature data.

[0129] The frequency characteristic synchronization data are processed to calculate the main frequency overlap and extract the cross-spectral coherence, resulting in main frequency overlap data and cross-spectral coherence data, respectively.

[0130] In one embodiment, the main frequency overlap is calculated as follows: ,in Main frequency overlap, The dominant frequency of the deformation trend signal. The dominant frequency of the microseismic signal. The main frequency average (for) );like This indicates "overlapping dominant frequencies"; cross-spectral coherence is calculated by using the Welch method to calculate the power spectral density of the two signals. Coherence is ,in Cross-spectral coherence (frequency) (place) The cross-power spectral density of the deformation signal and the microseismic signal. The power spectral density of the deformed signal, For the power spectral density of the microseismic signal, extract the maximum coherence value or the average coherence. ;like They believe that the spectral resonance synchronization is strong.

[0131] According to the main frequency overlap data and the cross spectrum coherence data, the coupling degree map of the deformation trend data and the microseismic monitoring data is constructed, and coupling degree map data is obtained.

[0132] In an embodiment, the system constructs a coupling degree function on a time sequence to measure the frequency domain coupling strength between the deformation trend signal and the microseismic signal at the same time. The coupling degree calculation adopts a weighted superposition strategy for definition, that is , wherein is the coupling degree function (changes with time), is the main frequency overlap weight coefficient, and the value is 0.5, is the main frequency overlap (a frequency domain similarity index), is the cross spectrum coherence weight coefficient, and the value is 0.5, is the cross spectrum average coherence. The system maps the above-mentioned time sequence coupling degree to a specific spatial measurement point position, binds the corresponding time sequence to the spatial coordinates according to the coordinate index of each measurement point, and forms a four-dimensional space-time coupling field tensor ; on the spatial structure, the system constructs a coupling graph network , each node is a spatial measurement point; the edge set is established according to the spatial adjacency relationship or the geological structure connection relationship between the measurement points; each node is additionally provided with a plurality of attribute fields, including the average coupling degree (time average in a certain period), the peak coupling degree (maximum value in a historical time period), the coupling frequency interval (main frequency statistical distribution result), and the signal cooperation duration.

[0133] The coupling degree map data is subjected to same frequency response extraction to obtain disturbance same frequency feature data.

[0134] In an embodiment, each node in the coupling degree map is traversed, and the following two indexes are used for judgment. First, the maximum coupling degree value of the node in the historical time sequence should be greater than 0.8, that is: , wherein represents the coupling degree score of the node at time , which is used to measure the synchronization degree between the microseismic frequency and the deformation trend frequency; second, the standard deviation of the main frequency of the node in a time window is less than 1.0 Hz, that is, the main frequency fluctuation is small, which shows the frequency stability: , the nodes satisfying the above two indicators will be marked as "co-frequency response points", indicating that there is strong frequency coupling and stable dominant frequency at this location within a certain period of time. Perform spatial aggregation operation on all nodes marked as co-frequency response points. Adopt three-dimensional voxel merging strategy to merge adjacent co-frequency response points with a distance not exceeding a certain threshold (such as 2 meters) into a continuous spatial region. Each merged region is taken as a response unit, and the following feature fields are extracted, including center coordinates, geometric center position of all points in the region; dominant frequency range, statistical result of dominant frequency value in the region (such as average dominant frequency and its interval); response duration, continuous high coupling time period from start to end (for example, more than 120 minutes above the threshold); regional coupling score, maximum or weighted average value of coupling degree in the region, to measure the resonance intensity of the region.

[0135] Optionally, the disturbance resonance identification comprises:

[0136] According to the disturbance co-frequency feature data, microseismic energy release density calculation and deformation gradient intensity calculation are performed, respectively obtaining microseismic energy release density data and deformation gradient intensity data;

[0137] In an embodiment, the microseismic data is voxelized to spatial unit grid as the center, the total energy of all microseismic events in the time window at this position is counted, and the total energy is normalized to the volume unit, wherein is the microseismic energy release density tensor, is the microseismic event index, is the voxel volume, is the release energy of the th event. According to the deformation trend tensor field , the spatial gradient is calculated: wherein is the spatial gradient of the deformation tensor field, is the partial derivative of the deformation tensor in the coordinate axis direction, is the partial derivative of the deformation tensor in the coordinate axis direction, is the partial derivative of the deformation tensor in the coordinate axis direction, and the modulus length is taken as the intensity index: wherein is the deformation gradient intensity index, is the spatial gradient of the deformation tensor field.

[0138] According to the microseismic energy release density data and the deformation gradient intensity data, light resonance judgment is performed to obtain resonance point data; ​

[0139] In an embodiment, a decision threshold is set, such as a microseismic energy density threshold , such as 100 J / m3; a deformation gradient intensity threshold , such as 0.05 mm / m; a same-frequency score threshold , output from the previous module, with a value such as 0.8. If a point satisfies the following conditions simultaneously: ， is the microseismic energy release density, is the energy threshold, with a value of 0.9, , is the deformation gradient intensity, is the deformation threshold, with a value of 0.05, and the same-frequency response score , is the same-frequency coupling score, is the coupling threshold, with a value of 200 J / m3, then the point is identified as a disturbance resonance point. Each resonance point is assigned the following information, including resonance energy weight , deformation response amplitude , start and end time, and spatial location.

[0140] According to the resonance point data, disturbance resonance data is generated to obtain the disturbance resonance data.

[0141] In an embodiment, based on spatial proximity, local density clustering (such as DBSCAN or spatial grid merging) is performed on the resonance point set to obtain a set of local resonance cores, each containing multiple resonance points. For each core region, the number of resonance points, the average microseismic energy density, the average deformation gradient, the main frequency interval aggregation degree, and the time concentration degree are counted to construct a multi-dimensional descriptor for describing the disturbance response characteristics of the region. The above information collectively constitutes a regional-level disturbance resonance descriptor for multi-dimensional description of the spatial distribution, energy characteristics, and time characteristics of the resonance aggregate.

[0142] Optionally, the resonance candidate unit extraction includes:

[0143] The disturbance resonance data is subjected to spatial continuity grid construction to obtain resonance point network data.

[0144] In an embodiment, the disturbance resonance point set of the previous module is input, each point having 3D spatial coordinates and resonance characteristic values such as energy density, resonance score, deformation gradient, etc. The underground space is divided into a voxel grid structure (such as each cube being a unit cell), and each voxel is numbered . All disturbance resonance points are assigned to the corresponding voxel unit : resonance point set, Construct a spatial topology graph based on 26 adjacency rules (up, down, left, right, front, back, and diagonal voxels). Each voxel is a graph node. Establish edges between adjacent voxels .

[0145] Local response filtering is performed on the resonance point network data to obtain response hotspot kernel data;

[0146] In one embodiment, the system uses each spatial voxel unit in the resonant point network structure. As a statistical unit, the following four key response indicators are extracted from the resonance point data, including the number of resonance points ( The total number of resonance points contained within the voxel; the average resonance intensity fraction ( : The arithmetic mean of the resonance scores of all resonance points in the voxel, used to measure the consistency and coupling strength of the disturbance in the region; average microseismic energy density ( ): Reflects the average intensity of energy release from local microseismic events; average deformation gradient intensity ( ): Represents the average level of local strain rate, revealing the severity of the structural response. To determine whether a voxel constitutes a response hotspot core, the system sets the following response screening conditions. If a voxel meets at least two of the following four indicators, it can be marked as a hotspot core voxel: Resonance point number threshold condition: (e.g., ≥3), for example, setting a threshold. The score is 3; Resonance intensity scoring criteria: For example, setting the threshold to 0.8; energy release intensity condition: For example, if the threshold is set to 100 J / m³; deformation gradient response conditions: For example, a threshold of 0.05 mm / m is set. For voxel elements that meet the above conditions, the system marks them as hotspot core voxels and outputs the following structured information as response hotspot core data, including the hotspot core number, assigning a unique identifier to each voxel that meets the conditions; spatial center coordinates, i.e., the geometric center position of the voxel element in three-dimensional space; and the number of resonance points (…). ); average resonance intensity fraction ( ).

[0147] Spatial topology filtering is performed based on the response hotspot kernel data to obtain resonance candidate unit data.

[0148] In one embodiment, the system is based on a three-dimensional voxel map structure. Perform topological connectivity analysis on each node in the graph. Each edge corresponds to a voxel unit containing a hotspot kernel. represents the spatial adjacency relationship between voxels (e.g., based on the 26-adjacency principle). The system identifies all connected subgraphs composed of hot spot voxels on the graph, each connected subgraph represents a set of voxels in the disturbance feature set that are mutually adjacent and have spatial continuity, and is regarded as a candidate resonance body unit. The system sets the following filtering rules: if the number of voxels contained in a connected subgraph is less than a preset threshold (e.g., 3 voxels), it is determined that the structure has weak stability and is rejected; the system sets the stability of the voxel neighborhood barycenter as the rejection criterion. The drift rate of the connected barycenter in the local window is calculated, and if the barycenter movement rate exceeds a set threshold (0.10-0.20), it is determined that the structure has strong disturbance noise and is not included in the effective candidate unit. For each connected resonance unit remaining, the system extracts the following aggregation attributes as its feature representation, including the spatial envelope volume, the minimum circumscribed envelope volume of the voxel set it occupies, which is used to quantify its spatial coverage; the highest score of the resonance point, which extracts the maximum value of the resonance intensity score in the contained voxels, which is used to measure the upper limit of the local disturbance intensity in the unit; the center coordinates (barycenter), which are determined by averaging the voxel center coordinates to determine the geometric barycenter of the resonance unit; the structure directionality and symmetry, which extract shape information such as the principal direction (e.g., the principal axis direction) and the distribution of symmetry axes.

[0149] Optionally, S3 includes:

[0150] S31, risk feature extraction is performed on the deformation disturbance resonance data to obtain risk feature data;

[0151] In an embodiment, the system receives deformation disturbance resonance data from a previous processing module, and each spatial unit in the data contains, for example, resonance intensity scores, which represent the coupling degree of deformation trend and microseismic frequency in the local area; deformation trend slope, which represents the deformation growth rate per unit time; microseismic energy release rate, which refers to the energy released by seismic events per unit time and unit volume. The system extracts the following four types of main feature indicators for each spatial unit, including maximum historical deformation, which reflects the extreme deformation amplitude during the historical monitoring period; current deformation rate, which represents the unit time change rate of deformation in the recent period; unit time microseismic event frequency, which counts the number of microseismic events in the current time window; unit volume microseismic energy density, which represents the energy release intensity of the local area; disturbance resonance intensity score, which combines the coincidence degree of the main frequency and the cross-spectral coherence (the weights are 0.5 and 0.5) to calculate the resonance performance score; resonance duration, which refers to the duration that meets the resonance condition; cross-spectral coherence between microseismic and deformation, which measures the frequency domain synchronization strength between the two types of signals.

[0152] S32, region risk classification is performed according to the risk feature data, to obtain region risk classification data;

[0153] In an embodiment, the system can use the following rule combination to make a hierarchical judgment on the space unit: if the resonance score , the current deformation rate , and the microseismic energy density , it is determined as high risk (Dangerous); if the resonance score and : it is determined as medium risk (Warning); if the above conditions are not met, it is defaulted as low risk (Normal); if some features are abnormal or missing, and the model confidence is lower than a set threshold, it can be marked as uncertain risk (Uncertain).

[0154] S33, level connectivity correction is performed according to the region risk classification data, to obtain risk connectivity graph data;

[0155] In an embodiment, the system performs grid division on the underground structure space, regards each grid unit as a node, and constructs a space adjacency graph , wherein represents a set of all grid nodes; represents a pair of nodes having a spatial adjacency relationship, i.e., two grids have a connection relationship of sharing edges, faces or corner points in three-dimensional space. In the adjacency graph structure, the adjacent nodes having the same risk level are aggregated to form several connected regions, which are called risk level connected clusters. Each connected cluster corresponds to a risk level region, has a unique identification number and a list of included grid units. For the extracted high-risk connected cluster, the system sets a minimum connected region judgment threshold, for example, a micro cluster with an area less than 10 square meters or containing less than 3 voxels. If a cluster meets the "small cluster" condition, the following processing logic is executed: the average risk level of the surrounding adjacent regions is calculated; if the adjacent regions are mainly medium or low risk, the high-risk region is determined as an isolated anomaly, and risk downgrading is performed; if the surrounding regions have high or very high risk levels, the original risk level is retained or weighted upgrading is performed to enhance connectivity. For the junction zone of multi-level risk regions, the system constructs an edge fusion region, and calculates the deformation gradient direction, deformation trend main axis and other indicators in the region to determine the dominant risk direction and risk level dominant attribute. According to the spatial gradient continuity or deformation trend consistency, the fusion region is assigned a unified risk level.

[0156] S34, stability layering is performed according to the risk connectivity graph data, to obtain risk region hierarchical data.

[0157] In an embodiment, the system builds a stability level standard based on the above multi-dimensional indicators, forming a hierarchical rule as follows (which can be defined as four levels): Level 1 (Level_1, extremely high risk): the region meets the following conditions at the same time: microseismic event is dense (microseismic density threshold, value is 8 times / m³), deformation is severe (deformation gradient threshold, value is 20 mm / m), disturbance intensity is high (disturbance intensity threshold, value is 0.7), and has large-span connected structure (connected span threshold, value is 15 m), which indicates the existence of instability zone; Level 2 (Level_2, higher risk): microseismic frequency is high (microseismic frequency threshold, value is 2 times / m³) or deformation rate is large (deformation rate threshold, value is 7 mm / m), but the connected range is small, and the stability is relatively good; Level 3 (Level_3, medium risk): microseismic events and deformation trends cross-distribute in the region (cross-distribution threshold, value is 2 times / m³), but the intensity is general; Level 4 (Level_4, low risk): disturbance and deformation are sparse (disturbance intensity threshold, value is 0.3), and the overall stability of the region is good. The above thresholds are set according to the accuracy of the monitoring system or expert experience. The above thresholds are set according to the accuracy of the monitoring system or expert experience.

[0158] Optionally, S4 includes:

[0159] S41, constructing a patrol candidate task unit graph according to the risk region classification data to obtain patrol candidate graph data;

[0160] ​​​​​​​​​​​​​​​​In an embodiment, the input data used in this step is the aforementioned risk zone classification data. In this data, each spatial unit (e.g. voxel unit or aggregated grid block) is attached with a risk level label, such as Level_1: very high risk; Level_2: high risk; Level_3: medium risk; Level_4: low risk. The system sets up inspection inclusion rules according to the risk level, and filters the areas that meet the conditions as candidate task unit nodes. The filtering rules are as follows: mandatory inclusion node, all spatial units marked as Level_1 (very high risk) and Level_2 (high risk) must be included in the inspection candidate graph as core inspection targets. Conditional inclusion node, for Level_3 (medium risk) areas, the system can set up a spatial buffer condition to determine whether to be included, and the setting rule is that if there is a Level_2 area within the preset spatial buffer radius (such as 10 meters), it will also be retained. Default exclusion node, Level_4 (low risk) areas are not included in the graph structure by default, unless they are the only passage unit to the high-risk area (for example, the only passage unit to the high-risk area), which can be retained as a connection edge. The system constructs a task unit graph with the filtered spatial units as the node set of the graph structure. Each node represents a spatial candidate task unit, which can be a grid block, a voxel aggregation area, or a scene logical task block. Edges are constructed according to the passage structure in the chamber design (such as the connection relationship of the roadway), or the unobstructed reachable path between two nodes in a three-dimensional Euclidean space is used as the connection standard. Each edge is assigned a weight value representing the passage cost. The weight can be defined as the pure geometric distance, or can be set as a composite weighted form that integrates risk attributes, such as: wherein is the edge weight between node and node , is the passage distance weighting coefficient, is the Euclidean distance between node and node , is the risk level weighting coefficient, is the risk level value of node , and is the risk level value of node .

[0161] S42, sorting the nodes of the inspection candidate graph data to obtain inspection sorted graph data;

[0162] In an embodiment, the system constructs a priority score function for each inspection candidate node as the basis for sorting. The priority score function is in the form of wherein is the priority score of node . is a risk level weight coefficient, taking a value of 0.6, is a risk level score of the node , is a position priority weight coefficient, taking a value of 0.2, is a position priority score of the node , is an access cost weight coefficient, taking a value of 0.2, is an access cost value of the node . The priority score function can set the weight coefficients by AHP (analytic hierarchy process) or experience parameter adjustment. According to the risk level and position segmentation (Level_1→Level_2→Level_3); in each risk level segment, the nodes are sorted according to the priority score; the nodes with lower access cost are preferentially selected in the same level segment.

[0163] S43, three-dimensional path planning construction is performed on the inspection sorting graph data, and three-dimensional path data is obtained;

[0164] In an embodiment, the system discretizes the cavern model by using a grid voxelization method in a three-dimensional space. The three-dimensional space is divided into equidistant cubic units (voxels), and each voxel point is taken as a position node in the graph; the graph structure is constructed according to the channel connectivity and physical constraints between adjacent voxels, and is used for path search. The path planning can select a three-dimensional A* algorithm, heuristic search, and efficient path finding combined with a heuristic function (such as Euclidean distance); a three-dimensional Dijkstra algorithm, shortest path search without heuristic conditions, to ensure a global optimal solution. The planning process takes the adjacent two nodes in the sorting task point as the path start point and end point, and performs path calculation pair by pair. In the path reachability judgment process, the system performs slope constraint to ensure the path executable in actual application. The system calculates the path slope between two adjacent voxel nodes and , and the formula is as follows: wherein is a path slope value, is an arctangent function (used for calculating the slope), is an elevation coordinate of the node , is an elevation coordinate of the node , is a horizontal coordinate (X-axis) of the node , is a horizontal coordinate (X-axis) of the node , is a horizontal coordinate (Y-axis) of the node , is a horizontal coordinate (Y-axis) of the node horizontal coordinate (Y axis) of the node , are the spatial coordinates of the nodes and respectively; if the slope value exceeds the critical angle set by the system (such as 45°), the connection is considered impassable and automatically excluded from the path graph. This rule can effectively exclude physically inaccessible paths due to excessive height difference or steep structure. The system performs smoothing processing on the original path node sequence. Methods such as B-spline or Catmull-Rom curve can be used. The three-dimensional path data output by the system is a structured path sequence, recording the spatial path planning results between all task points. Each path contains the starting task node number; the ending task node number; and a three-dimensional path point sequence (composed of consecutive coordinate points, representing the path trajectory).

[0165] S44, spatial constraint correction is performed on the three-dimensional path data to obtain inspection graph data for underground cavern inspection auxiliary operations.

[0166] In an embodiment, the system sets the following four types of space constraint types for the physical limitations and operational safety requirements in the underground environment where the inspection path is located, including if the path passes through a critical area of the structure (such as a dense area of surrounding rock cracks, a fault intersection area), path avoidance or replacement is required. Such areas can be determined by geological models such as structure risk maps, surrounding rock grade maps, etc.; if the path passes through an area where large equipment is arranged (such as ventilation equipment, supporting structure, storage facilities, etc.), the system identifies the equipment location through the facility shielding map or construction layout, and forces the path to avoid crossing and interfering; for the identified high-risk areas (such as microseismic active points, deformation mutation points, etc.), the system sets a minimum safety distance threshold, for example, 1.5 meters, requiring the inspection path to always maintain not less than the distance; in order to ensure the visibility and communication stability during the execution of the inspection task, the system requires the space where the path is located to be in the signal coverage area. The signal coverage information can be derived from the spatial distribution map of WiFi or ZigBee signal strength. The system checks each segment of the path data for space constraints, and if it finds that a segment of the path violates the above constraints, the following modification strategy is executed, if the space where the path segment is located still has available path space (such as open areas or lateral areas of gentle slopes), the system tries to find a feasible suboptimal path point near the original path and replace the path segment; if the path segment is completely in an impassable area (such as the back of a closed equipment, a signal blind area, or the inside of a high and steep slope), the path segment will be marked as "needs manual confirmation" to prompt the dispatcher to manually intervene or make on-site decisions. The modified inspection path data will be structured into three types of core layer content to form an inspection map, including a three-dimensional path point set, which records the point sequence of all inspection path segments (which can be a list of continuous spatial points or a combination of control points and smooth curves); path segment attributes, which label each path segment with attribute tags such as risk level (high / medium / low), traffic state (passage / avoidance / confirmation required), etc.; auxiliary layers output warning_points, which mark all positions that violate space constraints or manual confirmation points; safe_corridors, which identify and mark recommended path segments with good traffic continuity and low risk level. The output inspection map data can be spatially overlaid with a building information model (BIM) to form a visual inspection map, providing three-dimensional path guidance and risk avoidance assistance for actual operations.

[0167] Optionally, the application also provides an underground cavern inspection system for executing the underground cavern inspection method as described above, the underground cavern inspection system comprising:

[0168] a rock mass modeling and parameter inversion module, configured to obtain geological survey data, microseismic monitoring data, and historical deformation monitoring data; perform rock mass viscoelastic constitutive modeling according to the geological survey data and the microseismic monitoring data to obtain a rock mass viscoelastic model;

[0169] a deformation prediction and perturbation resonance analysis module, configured to perform deformation trend prediction according to a rock mass viscoelastic model and historical deformation monitoring data to obtain deformation trend data, and to perform deformation perturbation resonance analysis according to the deformation trend data and microseismic monitoring data to obtain deformation perturbation resonance data;

[0170] a multi-dimensional risk region grading module, configured to perform risk region grading on the deformation perturbation resonance data to obtain risk region grading data;

[0171] an intelligent inspection graph generation and task planning module, configured to perform inspection graph generation according to the risk region grading data to obtain inspection graph data for underground cavern inspection auxiliary operation.

[0172] Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting, the scope of the application being defined by the appended application file and not by the above description, so as to intend to encompass all variations falling within the meaning and scope of the equivalent requirements of the application file.

[0173] The above description is merely one specific implementation of the application, which enables those skilled in the art to understand or implement the application. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the application. Therefore, the application will not be limited to these embodiments shown herein, but will conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method of inspecting an underground cavern, characterized by, The method comprises: S1, acquiring geological survey data, microseismic monitoring data and historical deformation monitoring data; performing rock mass viscoelastic constitutive modeling according to the geological survey data and the microseismic monitoring data to obtain a rock mass viscoelastic model; S2, performing deformation trend prediction according to the rock mass viscoelastic model and the historical deformation monitoring data to obtain deformation trend data; performing deformation disturbance resonance analysis according to the deformation trend data and the microseismic monitoring data to obtain deformation disturbance resonance data; S3, classifying risk areas according to the deformation disturbance resonance data to obtain risk area classification data; S4, constructing a patrol candidate task unit graph according to the risk area classification data to obtain patrol candidate graph data; performing node sorting on the patrol candidate graph data to obtain patrol sorting graph data; performing three-dimensional path planning construction on the patrol sorting graph data to obtain three-dimensional path data; performing spatial constraint correction on the three-dimensional path data to obtain a patrol graph for assisting underground cavern patrol operations.

2. The method of claim 1, wherein, The S1 comprises: acquiring geological survey data, microseismic monitoring data and historical deformation monitoring data; performing constitutive model construction according to the geological survey data and the microseismic monitoring data to obtain a constitutive model; performing microseismic feature extraction according to the microseismic monitoring data to obtain microseismic feature data; performing viscous field inversion graph processing on the microseismic feature data to obtain viscous field inversion graph data; synchronously updating the constitutive model according to the viscous field inversion graph data to obtain a fracture network evolution model; performing coupling graph modeling on the fracture network evolution model and the geological survey data to obtain a rock mass viscoelastic model.

3. The method of claim 1, wherein, The deformation trend prediction comprises: performing deformation representation network structure construction according to the rock mass viscoelastic model and the historical deformation monitoring data to obtain deformation network representation data; performing trend change boundary identification on the deformation network representation data to obtain trend boundary data; performing variable response degree evaluation on the trend boundary data to obtain trend response degree data; performing trend mapping on the rock mass viscoelastic model according to the trend response degree data to obtain deformation trend data.

4. The method of claim 1, wherein, The deformation disturbance resonance analysis comprises: performing disturbance same-frequency feature extraction according to the deformation trend data and the microseismic monitoring data to obtain disturbance same-frequency feature data; performing disturbance resonance identification according to the disturbance same-frequency feature data to obtain disturbance resonance data; performing resonance candidate unit extraction on the disturbance resonance data to obtain resonance candidate unit data; performing clustering processing on the resonance candidate unit data to obtain resonance body data; performing resonance intensity classification on the resonance body data to obtain resonance classification data; performing mapping on the deformation trend data according to the resonance classification data to obtain deformation disturbance resonance data.

5. The method of claim 4, wherein, The disturbance same-frequency feature extraction comprises: performing empirical mode decomposition according to the deformation trend data and the microseismic monitoring data to obtain mode decomposition data; performing frequency feature extraction on the mode decomposition data to obtain frequency feature data; performing time synchronization on the frequency feature data to obtain frequency feature synchronization data; performing main frequency overlap degree calculation and cross-spectrum coherence degree extraction on the frequency feature synchronization data to obtain main frequency overlap degree data and cross-spectrum coherence degree data, respectively; According to the main frequency overlap data and the cross spectrum coherence data, coupling degree map data of the deformation trend data and the microseismic monitoring data are constructed. The coupling degree map data are subjected to same frequency response extraction to obtain disturbance same frequency characteristic data.

6. The method of claim 4, wherein, The disturbance resonance recognition includes: According to the disturbance same frequency characteristic data, microseismic energy release density calculation and deformation gradient intensity calculation are performed to obtain microseismic energy release density data and deformation gradient intensity data respectively. According to the microseismic energy release density data and the deformation gradient intensity data, light resonance judgment is performed to obtain resonance point data. According to the resonance point data, disturbance resonance data are generated to obtain disturbance resonance data.

7. The method of claim 4, wherein, The resonance candidate unit extraction includes: The disturbance resonance data are subjected to spatial continuity grid construction to obtain resonance point network data. The resonance point network data are subjected to local response screening to obtain response hot spot kernel data. According to the response hot spot kernel data, spatial topology filtering is performed to obtain resonance candidate unit data.

8. The method of claim 1, wherein, S3 includes: Risk feature data are obtained by extracting risk features from the deformation disturbance resonance data. Regional risk classification data are obtained by classifying regions according to the risk feature data. Risk connectivity graph data are obtained by correcting the grade connectivity according to the regional risk classification data. Risk area classification data are obtained by performing stability layering according to the risk connectivity graph data.

9. An underground cavern inspection system characterized by, The underground cavern inspection system for performing the underground cavern inspection method as claimed in claim 1 includes: The rock mass modeling and parameter inversion module is used to obtain geological survey data, microseismic monitoring data and historical deformation monitoring data; rock mass viscoelastic constitutive modeling is performed according to the geological survey data and the microseismic monitoring data to obtain a rock mass viscoelastic model. The deformation prediction and disturbance resonance analysis module is used to predict deformation trend according to the rock mass viscoelastic model and the historical deformation monitoring data to obtain deformation trend data; deformation disturbance resonance analysis is performed according to the deformation trend data and the microseismic monitoring data to obtain deformation disturbance resonance data. The multi-dimensional risk area classification module is used to classify risk areas according to the deformation disturbance resonance data to obtain risk area classification data. The intelligent inspection map generation and task planning module is used to generate an inspection map according to the risk area classification data to obtain inspection map data for underground cavern inspection auxiliary operation.

Images