Intelligent analysis method and system for underground structure based on multi-source geological data fusion

Abstract

This invention discloses an intelligent analysis method and system for underground structures based on multi-source geological data fusion, relating to the fields of engineering geology and geotechnical engineering. The method includes: S1, collecting raw geological data and applying a hierarchical clustering algorithm to group scale-scale fractures and interlayer slippage to obtain a preliminary stiffness distribution map of scale-scale rock mass units; S2, based on the preliminary stiffness distribution map of the scale-scale rock mass units, using finite element analysis to simulate the stress transfer path between units and determine the deformation mode and failure path of the scale-scale geological blocks. This intelligent analysis method and system for underground structures based on multi-source geological data fusion can accurately identify stress concentration, uneven settlement, and potential sliding surface locations in key areas such as weak interlayers in dam foundations, reducing judgment bias caused by empirical simplification. It provides more objective, stable, and repeatable analysis results for stability evaluation and risk prediction of underground engineering, and has good engineering applicability and promotional value.

Description

Technical Field

[0001] This invention relates to the field of engineering geology and geotechnical engineering technology, specifically to an intelligent analysis method and system for underground structures based on the fusion of multi-source geological data. Background Technology

[0002] Geological structure analysis is an important research area in earth science and engineering, with its findings widely applied in hydropower projects, transportation tunnels, underground storage facilities, and large-scale infrastructure construction. In these projects, the spatial distribution characteristics, mechanical properties, and deformation and failure behavior of underground rock mass structures under external loads directly affect the safety and long-term stability of the engineering structures. Therefore, how to rationally model and analyze underground structures has always been a key focus in the fields of engineering geology and geotechnical engineering.

[0003] In practical engineering, underground geological bodies are typically composed of multi-scale structures, including micro-fractures, inter-layer faulting, mesoscale rock mass units, and large-scale geological blocks. Significant mechanical coupling exists between these different scale structures. The distribution characteristics of small-scale structures affect the equivalent stiffness and strength of mesoscale rock mass units, while the combination of mesoscale units further determines the overall deformation pattern and potential failure path of the large-scale geological body. However, existing geological structure analysis methods often struggle to balance computational efficiency and the realism of mechanical responses when dealing with such multi-scale problems. On the one hand, some methods tend to retain as many detailed structures as possible, such as fractures and weak interlayers, during the modeling process to improve the geometric precision of the model. However, this results in a large model size and high computational complexity, making it unsuitable for regional or engineering-level numerical analysis. On the other hand, when simplifying complex geological bodies, existing methods often employ empirical or static simplification rules, ignoring the mechanical transmission relationships between different scale structures. This makes it difficult for the simplified model to accurately reflect the overall deformation characteristics and failure evolution of the original geological body under stress. Furthermore, in engineering analysis, the deformation patterns and potential failure locations of underground structures typically need to be verified by comparing numerical simulations with field monitoring data. However, existing technologies generally lack an analytical workflow that can establish an effective feedback mechanism between model construction, analysis, and verification. When there are significant deviations between simulation results and actual monitoring data, it is often necessary to rely on manual experience to adjust model parameters, lacking systematic iterative optimization methods, which makes it difficult to guarantee the stability and reliability of the analysis results. Especially in critical engineering parts such as dam foundations, issues such as weak interlayers, uneven settlement, and stress concentration have a significant impact on engineering safety. If the mechanical response characteristics of these critical parts cannot be reasonably represented in the simplified model in regional-scale analysis, it is easy to cause inaccurate judgment of the potential sliding surface location, thereby affecting the reliability of engineering safety assessment and risk prediction. Summary of the Invention

[0004] The purpose of this invention is to provide an intelligent analysis method and system for underground structures based on the fusion of multi-source geological data, thereby solving the problems existing in the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: an intelligent analysis method for underground structures based on multi-source geological data fusion, comprising: S1, collecting raw geological data and applying a hierarchical clustering algorithm to group scale fractures and interlayer faults to obtain a preliminary stiffness distribution map of scale rock mass units; S2, based on the preliminary stiffness distribution map of scale rock mass units, using finite element analysis to simulate the stress transfer path between units and determine the deformation mode and failure path of the scale geological block; S3, if the deformation mode of the scale geological block deviates from a preset threshold by more than a specified range, then the hierarchical clustering algorithm is reapplied by iteratively adjusting the grouping parameters to obtain an optimized stiffness distribution map of the scale rock mass units. S4. After obtaining the optimized scale rock mass unit stiffness distribution map, for key parts including the weak interlayer of the dam foundation, the finite element analysis method is used to calculate the local stress concentration and uneven settlement distribution to determine the location of the potential sliding surface; S5. When determining the location of the potential sliding surface, by integrating the deformation mode of the scale geological block with the local stress concentration data, a multi-scale nested simulation algorithm is applied to construct the overall mechanical response model to obtain the simplified prediction results of the macroscopic behavior of the underground structure; S6. The stiffness gradient and slip coordination relationship are extracted from the simplified prediction results of the macroscopic behavior of the geological structure, and the finite element analysis method is used to verify the matching degree with the actual monitoring data to determine the consistency of the model.

[0006] Preferably, step S1 includes obtaining raw geological data from field exploration, applying a hierarchical clustering algorithm to the geological data for scale-based fracture grouping, and obtaining fracture grouping results; performing inter-layer fault analysis based on the fracture grouping results and inter-layer fault data, and determining the fault-affected area by comparing displacement differences; integrating the fault-affected area with the fracture grouping results to divide the scale-based rock mass units and obtain unit boundary descriptions; obtaining stiffness parameters within the unit boundary descriptions, mapping preliminary stiffness distributions, and generating a draft distribution map; extracting stability indices from the draft distribution map, and optimizing the rock mass units to obtain a preliminary stiffness distribution map of the scale-based rock mass units.

[0007] Preferably, step S2 includes obtaining scale rock mass element parameters from a preliminary stiffness distribution map, generating a mesh structure using the finite element method to obtain an inter-element connection description; inputting geological load data to simulate stress transfer based on the inter-element connection description, and determining path distribution characteristics; calculating the block deformation response through the path distribution characteristics to determine the mode type; obtaining the strain threshold under the mode type, performing strain interaction analysis on the rock mass elements to obtain a failure path sequence; and integrating block stability indices based on the failure path sequence to determine the deformation mode and failure path of the scale geological block.

[0008] Preferably, step S3 includes obtaining threshold deviation detection results by comparing deformation modes with preset thresholds, performing iterative adjustment operations on deviation data through grouping parameters to obtain an adjusted parameter set; for the adjusted parameter set, using a hierarchical clustering algorithm to restart the grouping process of rock mass unit data to generate preliminary rock mass unit groups; based on the preliminary rock mass unit groups, performing an operation to obtain strain thresholds from geological load data to determine an inter-unit stress balance adjustment scheme; using the inter-unit stress balance adjustment scheme, performing simulated input processing on load data to determine changes in grid connection description; obtaining changes in grid connection description, integrating failure path sequences and stability indices through weighted calculations to obtain an optimized scale rock mass unit stiffness distribution map.

[0009] Preferably, step S4 includes obtaining material property parameters of the weak interlayer region of the dam foundation through an optimized scale rock mass element stiffness distribution map, thus obtaining local geological model data; constructing a mesh model using the finite element analysis method based on the local geological model data, and determining boundary conditions and load distribution; inputting hydraulic seepage influence parameters from the mesh model, which are pre-obtained from the geological exploration data of the dam foundation, and calculating the stress field changes within the weak interlayer of the dam foundation. The calculation process uses stress balance equations to integrate seepage pressure distribution and obtain local stress concentration values; integrating uneven settlement distribution data based on local stress concentration values, and adjusting model parameters if the settlement deviation exceeds a preset threshold to determine the settlement equilibrium state; obtaining the settlement equilibrium state and, combined with failure path simulation, determining the location of the potential sliding surface.

[0010] Preferably, step S5 includes obtaining local stress concentration data through the deformation mode of the geological block at different scales, integrating relevant parameters of the deformation mode, and determining the initial response distribution of the underground structure; for the initial response distribution, using a multi-scale nested simulation algorithm, constructing an overall mechanical response model by integrating relevant parameters of the deformation mode and local stress data, and obtaining a simplified macroscopic behavior distribution; obtaining relevant information about the geological block from the macroscopic behavior distribution, calculating the groundwater level influence parameters, and determining the preliminary location of the potential sliding surface; obtaining the preliminary location, adjusting the relevant parameters of the deformation mode in conjunction with the local stress concentration data, and determining the correction value for the sliding surface location; and integrating relevant data of the macroscopic behavior distribution based on the correction value to obtain the prediction result of the macroscopic behavior of the underground structure.

[0011] Preferably, step S6 includes: obtaining relevant data on stiffness gradient distribution from the geological structure prediction results; constructing a preliminary description of the mechanical response law based on the stiffness gradient distribution and the changing characteristics of underground deformation trends, thereby obtaining the initial distribution characteristics of the mechanical response law; obtaining the specific value range of the slip coordination parameter based on the initial distribution characteristics of the mechanical response law; adjusting the boundary conditions of the mechanical response law by mapping the slip coordination parameter with the local stress distribution, thereby determining the corrected response distribution range; integrating relevant information on macroscopic behavioral characteristics based on the corrected response distribution range; constructing an input dataset for the finite element analysis method under the influence of geological environmental variables, thereby obtaining preliminary comparison results of the finite element analysis method; obtaining the deviation value of the monitoring data comparison based on the preliminary comparison results; adjusting the weight ratio of the slip coordination parameter and the stiffness gradient distribution according to the data matching accuracy requirements, thereby determining the distribution range of the behavioral consistency value; if the distribution range of the behavioral consistency value meets a preset threshold, then combining the relevant data on local stress distribution and underground deformation trends to generate the final geological structure prediction adjustment scheme and determine the verification conclusion of the macroscopic behavioral characteristics.

[0012] Preferably, the method further includes S7: For the determined model coordination consistency, if the matching degree is lower than a preset threshold, the scale fractures and interlayer faults are regrouped through a feedback loop to obtain the final regional geological evaluation model. Specifically, this includes obtaining specific data on fault reconstruction correlation from the fracture grouping iteration. The fracture grouping iteration performs multiple iterative classifications using the geometric distribution and stress field data of scale fractures. The fault reconstruction correlation calculates the correlation strength using the interlayer displacement vector and friction coefficient. Based on the results of the matching degree threshold verification, an initial framework for feedback loop optimization is constructed to obtain the preliminary distribution of model consistency adjustment. For the preliminary distribution, relevant information from scale feature simulation is integrated. The scale feature simulation uses the finite difference method to process fracture density and permeability data. Through the process of interlayer interaction correction, the interaction vector is adjusted according to the displacement gradient and shear modulus. The boundary conditions of the deviation value correction process are adjusted to determine the correction range of the integrated geological evaluation.

[0013] Preferably, step S7 further includes obtaining the input dataset for constructing the regional model based on the corrected value range, wherein the input dataset contains the corrected stress distribution and deformation modulus, and under the influence of the coordination parameter fusion, wherein the coordination parameter fusion combines the stiffness gradient and slip parameter through a weighted average method, performs the deformation trend calibration operation, and obtains the comparison result of the final scheme optimization; through the comparison result, the range of the behavior consistency value is determined, and if it is lower than the preset threshold, the fracture grouping iteration and fault reconstruction association are re-executed, wherein the re-execution includes updating the stress field data and friction coefficient, and obtaining the final regional geological evaluation model.

[0014] A multi-source geological data fusion-based intelligent analysis system for underground structures is used to implement the steps of the aforementioned multi-source geological data fusion-based intelligent analysis method for underground structures. The system includes a data acquisition and grouping module, which acquires raw geological data and applies a hierarchical clustering algorithm to group scale-scale fractures and inter-layer faults to obtain a preliminary stiffness distribution map of scale-scale rock mass units; a deformation mode analysis module, which, based on the preliminary stiffness distribution map of scale-scale rock mass units, uses finite element analysis to simulate the stress transfer path between units and determine the deformation mode and failure path of the scale-scale geological block; an iterative optimization module, which, if the deformation mode of the scale-scale geological block deviates from a preset threshold by more than a specified range, iteratively adjusts the grouping parameters and reapplies the hierarchical clustering algorithm to obtain an optimized stiffness distribution map of the scale-scale rock mass units; and a local mechanical response analysis module, which obtains the optimized stiffness distribution map of the scale-scale rock mass units. After mapping the stress distribution, for key areas including weak interlayers in the dam foundation, finite element analysis is used to calculate local stress concentration and uneven settlement distribution to determine the location of potential sliding surfaces. The overall mechanical response modeling module, when determining the location of potential sliding surfaces, integrates the deformation patterns of scale geological blocks with local stress concentration data, and applies a multi-scale nested simulation algorithm to construct an overall mechanical response model, obtaining simplified predictions of the macroscopic behavior of the underground structure. The consistency verification module extracts stiffness gradients and slip coordination relationships from the simplified macroscopic behavior predictions of the geological structure, and uses finite element analysis to verify the matching degree with actual monitoring data, determining the model's consistency. The feedback evaluation module, for the determined model consistency, if the matching degree is lower than a preset threshold, regroups scale fractures and interlayer faults through a feedback loop to obtain the final regional geological evaluation model.

[0015] As can be seen from the above technical solution, the present invention has the following beneficial effects: This intelligent analysis method and system for underground structures, based on the fusion of multi-source geological data, can balance computational efficiency and the realism of mechanical response in regional-scale underground structure analysis. By grouping fractures and interlayer slippage at different scales and combining finite element analysis with multi-scale nested simulation, it realizes the mechanical transfer path analysis from mesoscale rock mass units to large-scale geological blocks, effectively reflecting the mechanical coupling relationship between multi-scale geological structures. Simultaneously, by comparing and verifying deformation patterns and stress concentration results with real monitoring data and introducing a feedback iteration mechanism to dynamically correct the model, it improves the consistency and reliability between underground structure analysis results and actual engineering behavior. This method can accurately identify stress concentration, uneven settlement, and potential sliding surface locations in key areas such as weak interlayers in dam foundations, while simplifying the model. It reduces judgment bias caused by empirical simplification, providing more objective, stable, and repeatable analysis results for underground engineering stability evaluation and risk prediction, and has good engineering applicability and promotional value. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the overall process of the intelligent analysis method for underground structures according to the present invention. Figure 2 This is a block diagram of the system structure of the present invention. Detailed Implementation

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

[0018] like Figure 1 As shown, this invention provides a technical solution: an intelligent analysis method for underground structures based on multi-source geological data fusion, comprising: S1. By collecting raw geological data and applying a hierarchical clustering algorithm to group scale fractures and interlayer faults, a preliminary stiffness distribution map of scale rock mass units is obtained. S2. Based on the preliminary stiffness distribution map of the scale rock mass unit, the finite element analysis method is used to simulate the stress transfer path between units and determine the deformation mode and failure path of the scale geological block. S3. If the deformation pattern of the scale geological block deviates from the preset threshold by more than a specified range, the hierarchical clustering algorithm is reapplied by iteratively adjusting the grouping parameters to obtain an optimized scale rock mass unit stiffness distribution map. S4. After obtaining the optimized scale rock mass unit stiffness distribution map, for key parts including the weak interlayer of the dam foundation, the finite element analysis method is used to calculate the local stress concentration and uneven settlement distribution, and to determine the location of potential sliding surfaces. S5. When determining the location of potential sliding surfaces, by integrating the deformation patterns of scale geological blocks and local stress concentration data, a multi-scale nested simulation algorithm is applied to construct an overall mechanical response model, and a simplified prediction result of the macroscopic behavior of underground structures is obtained. S6. Extract stiffness gradient and slip compatibility from the simplified macroscopic behavior prediction results of geological structures, and use finite element analysis to verify the matching degree with the actual monitoring data to determine the consistency of the model. S7. For the determined model coordination consistency, if the matching degree is lower than the preset threshold, the scale fractures and interlayer faults are regrouped through feedback loop to obtain the final regional geological evaluation model.

[0019] In the above implementation method, this method is based on multi-source geological data, including borehole data, geological logging information, seismic exploration data, and in-situ monitoring data. By uniformly modeling geological information from different sources and scales, a comprehensive analysis of underground structures is achieved. First, a hierarchical clustering algorithm is used to group scale-specific fractures and interlayer faults, merging structural surfaces with similar mechanical properties into rock mass units of the same scale, thus forming a preliminary stiffness distribution description on a macroscopic level. Subsequently, based on this stiffness distribution map, the finite element analysis method is introduced to numerically simulate the stress transfer paths between scale-specific rock mass units, revealing the overall deformation mode and potential failure paths of the geological block under external loads.

[0020] When the simulated deformation pattern deviates significantly from the preset threshold, the grouping parameters in the hierarchical clustering algorithm are iteratively adjusted to make the rock mass unit division more consistent with actual geological conditions, thereby continuously correcting the stiffness distribution map. Based on this, for key components such as weak interlayers in the dam foundation, a locally refined finite element model is used to calculate stress concentration and uneven settlement distribution, and the location of potential sliding surfaces is determined in conjunction with rock mass structural characteristics. Furthermore, a multi-scale nested simulation algorithm is used to integrate the deformation patterns of scale geological blocks with the results of local stress concentration, constructing a unified overall mechanical response model to predict the macroscopic behavior of underground structures. Finally, the predicted results are compared and verified with actual monitoring data to evaluate the consistency of the model in terms of stiffness gradient and slip coordination. The model parameters are continuously corrected through feedback loops until a stable and reliable regional geological evaluation model is obtained.

[0021] In the above implementation, by introducing the coupled application of hierarchical clustering algorithm and finite element analysis method, the effective fusion of multi-scale and multi-source geological data is achieved, enabling a more realistic characterization of the stiffness distribution and mechanical behavior of underground structures. This method significantly improves the model's adaptability to actual geological conditions through an iterative optimization mechanism, helping to reduce errors caused by unreasonable artificial partitioning in traditional analysis methods. Simultaneously, focused analysis on key components such as weak interlayers in the dam foundation leads to more accurate identification of potential sliding surfaces, thereby enhancing the reliability of underground structure stability evaluation. Verification through matching with real monitoring data further strengthens the credibility of the analysis results, providing strong technical support for the safety assessment and risk control of major projects.

[0022] S1 includes obtaining raw geological data from field exploration, applying a hierarchical clustering algorithm to group the geological data into scale fractures, and obtaining fracture grouping results; performing inter-layer fault analysis based on the fracture grouping results and comparing displacement differences to determine the fault-affected areas; integrating the fault-affected areas with the fracture grouping results to divide the scale rock mass units and obtain unit boundary descriptions; obtaining stiffness parameters within the unit boundary descriptions, mapping preliminary stiffness distributions, and generating a draft distribution map; extracting stability indices from the draft distribution map, optimizing the rock mass units into groups, and obtaining a preliminary stiffness distribution map of the scale rock mass units.

[0023] In this implementation, the first step is to collect and standardize the raw geological data. The raw geological data consists of borehole information, core logging records, surface and cavern exposure mapping results, geophysical interpretation results, in-situ testing results, and displacement monitoring results. To address the issues of inconsistent coordinate benchmarks, varying sampling densities, and differences in recording granularity among data from different sources, a unified spatial reference system was first established. Specifically, the following steps were taken: First, the coordinate and elevation benchmarks used in the project were determined, and all borehole openings, survey lines, measuring points, and monitoring points were converted to three-dimensional coordinates under the same benchmark. Second, for cases where multiple records exist at the same location, they were prioritized according to chronological order and data reliability, retaining the highest-priority record and establishing a traceable replacement relationship. Third, significant anomalies were identified based on source consistency and spatial continuity, characterized by sudden jumps in values ​​within the same geological body that are significantly inconsistent with the surrounding area, or sudden changes in values ​​at the same measuring point that are incompatible with instrument accuracy within a short period. After removing anomalies, the missing data was supplemented using spatial interpolation results from neighboring measuring points. Fourth, missing segments were filled in, based on statistically representative values ​​of the same stratum and lithology, and verified through comparison with adjacent boreholes to ensure that the supplemented results did not introduce new illusions of dense structural planes or weak zones. Key parameters involved in the above unification process included: data priority, anomaly identification threshold, and missing data filling range. Data priority is determined by the verifiability of the data acquisition method, with in-situ testing and direct exposure records having higher priority than indirect interpretation results. Anomaly identification thresholds are jointly determined by instrument accuracy, historical stable fluctuation range, and spatial continuity constraints. First, the upper bound of error is given by the calibration accuracy of the monitoring system, and then checked by the upper bound of fluctuation statistics under stable operating conditions. The final threshold is the more stringent of the two. The missing data completion range is determined by the maximum distance between adjacent valid data points. If the distance exceeds this range, supplementary exploration will be triggered or the confidence level of the area will be reduced, thereby avoiding over-inference.

[0024] After data unification, the calculation process for scale-based fracture grouping begins. Fracture grouping is based on the structural similarity of fracture sets. Fracture similarity is jointly characterized by fracture spatial location features, orientation features, spacing density features, extension scale features, continuity features, and filling and weathering features. To ensure comparability of features with different dimensions in the same clustering calculation, each feature is first normalized using linear scaling across the entire statistical area, ensuring that each feature participates in distance calculation within a uniform scale, thus preventing any single feature from dominating the grouping results due to excessive size. Subsequently, a pairwise fracture similarity matrix is ​​constructed. The similarity calculation follows the rule that "the smaller the azimuth difference, the smaller the spacing difference, the smaller the extension scale difference, and the higher the spatial proximity, the higher the similarity." For azimuth features, periodically consistent angle differences are used for calculation to maintain continuity when the orientation crosses boundaries. For spatial proximity, three-dimensional distance is used and stratigraphic interface constraints are introduced to prevent fractures across strata from being misclassified as belonging to the same group due to similar projections. The hierarchical clustering merging strategy adopts the principle of "minimizing the increase in intra-group differences" for merging step by step. During the step-by-step merging process, the incremental curve of intra-group differences brought about by each merging is recorded. When the incremental curve shows a significant leap, it indicates that continuing to merge will forcibly merge fractures with significant differences in structural characteristics. The grouping is truncated at this level, and the fracture grouping results are output. The key parameters involved in this stage include: feature set, feature weight, similarity calculation rules, and truncation level threshold. The feature set is determined by data availability and engineering sensitivity. If the project is significantly controlled by structural surfaces, the weight of orientation and continuity features is increased. Feature weights are determined by two types of criteria: the first type is the statistics of instability control factors in historical engineering areas, and the second type is the correlation analysis between the monitoring response and structural characteristics in this area. The stronger the correlation, the higher the weight. The truncation level threshold is determined by the leap point of the intra-group difference incremental curve and is verified using the spatial continuity of the fracture group. The verification criterion is that the fractures in the same group form a continuous or quasi-continuous banded distribution in space. If a large-scale discrete point distribution occurs, it is regressed to the previous level.

[0025] After obtaining the fracture grouping results, the process of inter-layer fault analysis and fault influence zone determination begins. Fault analysis uses inter-layer relative displacement as the core quantity, calculated from the displacement difference between upper and lower layers and between structural blocks on both sides within the same profile or monitoring sequence. Specifically, time synchronization and noise suppression are first performed on the displacement data. Time synchronization is achieved through a unified sampling period and interpolation alignment, while noise suppression is achieved through a sliding window smoothing process. The window length is determined by both the monitoring sampling frequency and the frequency of engineering loading changes, ensuring that the true deformation trend is preserved without excessive smoothing. Subsequently, the inter-layer displacement difference is calculated at each spatial location, and the rate of change of the displacement difference is calculated along the spatial direction to identify the start and end boundaries of the fault zone. When the displacement difference exceeds the fault identification threshold and forms a continuous area in space, the area is determined to be a fault influence zone. The determination of continuous areas is achieved through connectivity analysis. Connectivity analysis requires that the displacement difference between adjacent grids or adjacent measuring points all exceed the threshold and the continuous length exceeds the minimum continuous scale. The minimum continuous scale is determined by the structural dimensions of key engineering components and the spacing between monitoring points, ensuring that the identification results correspond to real geological zones rather than isolated noise points. The determination of the fault identification threshold adopts a "stable period statistical upper bound" strategy. First, a period in which the project is in a stable state is selected as the baseline period. Within this baseline period, the mean and dispersion of inter-story displacement differences are statistically analyzed. The threshold is calculated by superimposing the mean of the baseline period with a high-confidence upper bound for dispersion, and then verifying it against historical extreme values ​​of the baseline period. The final threshold is not lower than the upper bound of historical extreme values, thus ensuring that fault identification targets abnormal deformation rather than normal fluctuations. Key parameters involved in this stage include: smoothing window length, fault identification threshold, and minimum continuous scale of continuous areas. The window length is determined by the sampling frequency and loading change period; the fault identification threshold is determined by the stable period statistical upper bound and verified against historical extreme values; the minimum continuous scale is jointly determined by the influence range of key structures and the monitoring point density.

[0026] After determining the fault-affected area, the process involves integrating the fault-affected area with the fracture grouping results to delineate rock mass units at different scales. The core of this integration lies in using the structural control boundary formed by the fracture grouping as the initial partitioning framework, and the deformation control boundary formed by the fault-affected area as a strong constraint correction term. The two are then superimposed to form the final unit boundary description. Specifically, the fracture grouping results are first projected onto a unified spatial grid or polygonal partitioning framework to form fracture control partition surfaces. Then, the fault-affected area is transformed into an influence zone surface within the same framework, and its overlap with the fracture control partition surface is calculated. When the fault influence zone crosses a fracture partition and the overlap ratio exceeds the overlap threshold, partition correction is triggered. The correction method involves dividing the fracture partition along the boundary line of the fault influence zone, ensuring that the boundaries of the sub-units after division are consistent with the fault zone boundaries, thereby incorporating potential slip and discontinuous deformation into the unit boundary expression. The overlap threshold is determined by the engineering sensitivity to fault displacement. This sensitivity is determined by the correspondence between the monitoring response and the location of the fault zone. If the monitoring response exhibits abrupt changes in displacement or settlement near the fault zone, a more stringent overlap threshold is applied to strengthen the control of the fault zone over the element boundaries. After segmentation, element boundary descriptions are output. These descriptions include the spatial extent of each element, its adjacency relationships, and its contact relationship with the fault zone, providing a clear index for subsequent parameter extraction.

[0027] The process then proceeds to acquire stiffness parameters and map them to form a draft distribution map. Stiffness parameters are based on quantities reflecting the rock mass's resistance to deformation, and their sources include indoor mechanical test results, in-situ loading test results, wave velocity-based conversion results, and inversion-based mechanical parameter identification results. To ensure that parameters from multiple sources within the same unit can be fused at the numerical level, a consistency check is first performed on each source parameter. The consistency check is based on the correspondence between lithological type, degree of structural plane development, water-bearing state, and weathering degree. If a parameter is significantly inconsistent with its corresponding geological attribute, a source review or weight reduction is triggered. Subsequently, parameter aggregation is performed on rock mass units at each scale. During aggregation, representative stiffness values ​​of the units are calculated according to weights, which are determined by data reliability, spatial representativeness, and sampling coverage. Data reliability is determined by the acquisition method, with direct testing being higher than indirect conversion, and indirect conversion being higher than purely empirical values. Spatial representativeness is determined by the distance between the test point and the geometric center of the unit, as well as its relative position to the main structural surfaces; the closer the test point is to the unit core and the more it is located outside the main structural surfaces, the higher its representativeness. Sampling coverage is determined by the number of effective test points within the unit and the area coverage ratio; the higher the coverage ratio, the higher the weight. If there are obvious zoning characteristics within a unit, manifested as stable high and low stiffness partitions within the same unit, a secondary mapping within the unit is triggered. Representative values ​​of the partitions are used to assign zoning values ​​to the unit, so that the draft distribution map reflects the main stiffness variation trend within the unit. After completing the unit-level assignment, the unit stiffness values ​​are mapped to the entire spatial framework to form a draft stiffness distribution map.

[0028] After drafting the stiffness distribution map, the process of extracting and optimizing stability indices begins, aiming to elevate "reasonable geometric zoning" to "coordinated mechanical expression." The calculation of stability indices revolves around three types of issues: abrupt stiffness jumps between adjacent elements, uneven stiffness within elements, and incoordination between element boundaries and slip zones. Abrupt stiffness jumps between adjacent elements are obtained by calculating the difference in representative stiffness values ​​between elements on either side of a shared boundary, and weighted by the boundary length, making jumps at longer boundaries more significant in the evaluation. Uneven stiffness within elements is obtained by statistically analyzing the dispersion of mapping values ​​within the element; a dispersion exceeding a threshold indicates the existence of zonal structures within the element that are not represented by the boundary. Incoordination between boundaries and slip zones is obtained by calculating the deviation between the element boundary and the boundary of the slip zone; a deviation exceeding a threshold indicates that the element boundary does not follow the main discontinuous deformation zone. The thresholds mentioned above are all determined using a "regional statistical boundary" strategy: First, statistical distributions are formed for the jump values, discrete values, and deviation values ​​across the entire region, identifying the boundary points between concentrated intervals and anomaly tails, and using these boundary points as initial threshold values. Then, historical risk points in key engineering locations are used for verification. If a historical risk point is not identified by the threshold, the threshold is adjusted towards a stricter direction until all risk points are covered. After completing the index calculation and threshold determination, group optimization is performed: for regions where the stiffness difference between adjacent units has long been concentrated and the slip zone constraint is not significant, merging is performed to match the unit scale with mechanical consistency; for regions where the internal dispersion of a unit exceeds the threshold or the boundary deviates from the slip zone boundary by more than the threshold, segmentation is performed to ensure that the unit boundary responds to the main controlling discontinuous deformation and the main controlling stiffness zoning. After optimization, the stability index is recalculated. If all indices fall within the concentrated interval and the key location verification passes, a preliminary stiffness distribution map of the scale rock mass unit is output.

[0029] S2 includes obtaining scale rock mass element parameters from the preliminary stiffness distribution map, generating a mesh structure using the finite element method, and obtaining an inter-element connection description; for the inter-element connection description, inputting geological load data to simulate stress transfer and determine path distribution characteristics; calculating the block deformation response through path distribution characteristics to determine the mode type; obtaining the strain threshold under the mode type, performing strain interaction analysis on the rock mass elements to obtain the failure path sequence; and integrating the block stability index based on the failure path sequence to determine the deformation mode and failure path of the scale geological block.

[0030] In this embodiment, when obtaining scale rock mass element parameters from the preliminary stiffness distribution map, the geometric range of the element (directly determined by the three-dimensional closed boundary defined by the element boundary description) and the element volume (obtained by spatial integration of the three-dimensional boundary, and volume consistency verification is performed at the shared boundary of adjacent elements to avoid overlap or voids) are first determined based on the element boundary description. Then, the representative stiffness of the element (obtained by aggregating the values ​​assigned in the preliminary stiffness distribution map within the spatial range of the element, using a volume-weighted average method to ensure that sub-regions with high volume proportions contribute more to the representative value) and the element stiffness zoning are determined (obtained by identifying stiffness abrupt change zones within the same element in the preliminary stiffness distribution map; the threshold for abrupt change zone identification is determined simultaneously during identification, and the threshold is taken as the abnormal boundary point of the stiffness difference distribution in the entire area, and verification is performed at key locations to ensure that abrupt change zones correspond to real structural control rather than discrete noise). Simultaneously, the unit density (determined by indoor density testing and lithological data; if multiple values ​​exist, the weight is determined by sampling coverage, which is determined by the ratio of the number of sampling points within the unit to the unit volume) and the unit Poisson correlation value (determined by indoor mechanical tests or in-situ loading tests; if the test values ​​are discrete, they are screened based on lithological consistency and loading path consistency, with priority given to units of the same lithology, same water-bearing state, and same stress level) are extracted. For units with directional characteristics, directional parameters are determined (obtained from the statistical analysis of the orientation of dominant structural planes; dominant structural planes are determined by the group with the highest spatial continuity and the highest proportion of numbers among the fracture groups; isolated fracture points are removed during the statistical analysis to avoid the directionality being influenced by anomalies). For unit boundaries with weak interlayers or structural surfaces, determine the interface mechanical parameters (determined by the roughness of the structural surface, the properties of the filling material, the water content, and historical shear test data; friction parameters are taken from the median value of the test and reduced according to the water content, the reduction range is determined by the sensitivity evaluation of the water content to shear strength, and the sensitivity evaluation is determined based on the difference between dry and wet comparison tests of similar materials; bonding parameters are jointly determined by the strength grade of the filling material and the opening of the structural surface, the opening is obtained from geological logging statistics and the unfavorable quantile value is taken to ensure the expression of the safe side of the interface).

[0031] After extracting the element parameters, the finite element method is used to generate the mesh structure and form the inter-element connection description. Before generating the mesh, the outer boundary of the computational domain is determined (determined by the engineering analysis scope, covering key parts and their influence zones; the boundary of the influence zone is jointly determined by the load propagation distance and the influence range of the monitored anomalies; the propagation distance is obtained by statistically analyzing the high-stiffness channel length of the preliminary stiffness field, and the influence range of the monitored anomalies is obtained by statistically analyzing the spatial expansion boundary of historical anomalies). Then, mesh generation is performed: a volume mesh is established for each scale of rock mass element. The volume mesh size is determined by the minimum structural characteristic scale and the load gradient change scale; the minimum structural characteristic scale is taken as the minimum value of the weak interlayer thickness, the spacing between the main structural surfaces, and the width of the fault influence zone; the load gradient change scale is obtained by statistically analyzing the spatial variation amplitude of the initial stress field and the external load additional stress; a more stringent scale is taken to ensure that key geometric and stress changes are distinguished. Local densification is performed on key components. The scope of key components (determined by the location of weak interlayers in the dam foundation, fault influence zones, and high-risk areas of stress concentration, with the location of high-risk areas obtained by superimposing historical monitoring anomalies and geologically weak zones) and densification ratio (determined jointly by the upper limit of the computational scale and convergence stability requirements; convergence stability is judged by the change in results before and after mesh densification, with the change threshold determined simultaneously during the judgment; the threshold is the boundary value at which the difference in key response quantities under two subdivisions enters the stable interval). To form a description of inter-unit connections, interface consistency processing is performed on the shared boundaries of adjacent units: if a continuous force transmission boundary is used, nodes on the shared boundary are aligned (one-to-one corresponding node pairs are generated by boundary geometric projection mapping and error correction is performed; the upper limit of error is determined by the mesh size, and exceeding the upper limit triggers boundary re-subdivision); if a contact force transmission boundary is used, matching contact patches are formed on the shared boundary (generated by boundary triangulation, with the area difference between the two patches controlled within a quality threshold; the quality threshold is the anomalous boundary point of the patch area distribution and is locally corrected). The connection description includes the relationship between adjacent element numbers (uniquely determined by the shared boundary index), interface area (obtained by summing the area of ​​patches), interface normal direction (determined by the boundary geometry external normal consistency rule, which uses the external normal as a unified reference to avoid direction reversal), interface constraint type (determined by whether the interface corresponds to a structural surface or weak interlayer, based on the geological elements crossed by the boundary), and interface frictional bonding properties (directly assigned by the aforementioned interface mechanical parameters and adjusted in conjunction with the water content and filling type, with the adjustment rule prioritizing unfavorable combinations to avoid overestimating shear resistance). Mesh quality control is implemented throughout the meshing process. Quality indicators include minimum interior angle (used to avoid numerical ill-conditioning caused by thin, sharp elements; its threshold is obtained statistically from the numerical solution stability interval), element aspect ratio (used to limit overstretched elements; its threshold is obtained by assessing the sensitivity of key component responses to aspect ratio), and element distortion degree (used to limit geometric distortion; its threshold is determined by the boundary point where the element Jacobian discriminant enters the stability interval). If any indicator falls below the threshold, local re-meshment is triggered until the constraints are satisfied.

[0032] After obtaining the inter-unit connectivity description, geological load data is input and stress transfer simulation is performed. The geological load data is first processed by: self-weight load (determined by unit density and gravity direction, with the gravity direction defined by the coordinate datum and the density taken as the representative value of the aforementioned unit density) is applied to the entire volume unit; initial geostress (determined by regional stress testing, inversion identification, or existing engineering data, with values ​​selected using a multi-source consistency screening method, the screening criteria being that directionality and magnitude remain continuous within similar geological bodies, and in the event of abrupt changes, the source with higher reliability is prioritized, and the difference area is recorded as an uncertainty area) is applied through the initial stress field method; water pressure and seepage pressure (determined by groundwater level and pore pressure monitoring, with water level determined by long-term statistical high water level control values, and pore pressure determined by representative values ​​from the stable segment of the monitoring point time series, and values ​​assigned to seepage channels and low-permeability zones to maintain spatial gradient consistency with the geological structure) are applied. The loads are applied to the corresponding boundaries and volume elements; unloading and reloading caused by the construction stage (determined by aligning construction records and monitoring response inflection points, with stage division boundaries determined by the time point when the monitoring curve shows a trend reversal, and the reversal threshold being the upper limit of the stable segment fluctuation plus a safety margin) are applied gradually in stages; the volume deformation effect caused by temperature (determined by statistical analysis of temperature field measurement points, taking the unfavorable temperature difference combination within the design period and converting it according to the material thermal response parameters, with the conversion parameters determined by material tests or empirical data and verified with the monitoring thermal displacement response) are input synchronously; the inertial effect caused by seismic action (determined by the design seismic conditions, applied using representative time histories or equivalent inertial fields, with representativeness selected based on the matching degree between the site category and the design response spectrum, the matching degree being based on the energy coverage of the key frequency band and selecting unfavorable combinations). The load combination sequence follows the process of first establishing the initial stress field and then superimposing external loads, with the external load superposition sequence consistent with the construction stage, ensuring that the model evolution is consistent with the actual loading path. Boundary constraints are determined in conjunction with the location of the outer boundary: the strength of the far-field boundary constraint is determined by the relative scale of the distance from the outer boundary to the key part, the scale is the ratio of the distance from the outer boundary to the characteristic scale of the key part, and is verified by boundary sensitivity. The verification is judged by the change in the key response before and after boundary adjustment, and the threshold of the change is the dividing point of entering the stable interval.

[0033] After the stress transfer simulation is completed, the path distribution characteristics are extracted and a traceable transfer channel description is formed. The extraction of path distribution characteristics first identifies the dominant high-value zone from the global stress field: the dominant high-value zone (determined by the continuous high-value region formed by the stress field in space, the high-value judgment threshold is taken as the high quantile value of the stress distribution and checked at key parts, the verification basis is whether the high-value zone crosses the known weak zone and the monitoring anomaly area, if not covered, the threshold is adjusted to a more stringent direction until it is covered) is tracked along the continuous direction in space. During the tracking, cross-unit continuity is required (determined by the continuity of force transmission intensity on the shared boundary of adjacent units, the continuity judgment threshold is taken as the upper limit of the concentrated interval of the interface force transmission intensity distribution). Then the interface force transmission intensity is calculated (characterized by the normal contact pressure, tangential force transmission level and slip trend on the shared boundary, the normal contact pressure and tangential force transmission level are obtained by integrating the reaction force distribution on the interface by finite element solution, the slip trend is obtained by the cumulative amount of relative displacement on both sides of the interface in the tangential direction, and confirmed by the monotonicity of the time or stage sequence, the confirmation threshold is taken as the upper limit of the monitoring noise to exclude numerical oscillation). The concentration index for local concentration effects is calculated (determined by the abrupt increase in the local peak value relative to the representative value in the neighborhood, the neighborhood range being jointly determined by the grid size and structural feature scale; the abrupt increase threshold is taken as the abnormal boundary point of the abrupt increase distribution across the entire domain, and verified by grid refinement; the verification standard is that the concentration positions remain consistent before and after refinement, and the intensity ranking remains stable). The transmission attenuation across weak interlayers or structural surfaces is calculated (obtained by comparing the channel intensity on both sides of the interface; the comparison range is determined by a fixed distance band on both sides of the interface normal; the thickness of the distance band is determined by the grid size; the attenuation judgment threshold is taken as the abnormal boundary point of the attenuation distribution and verified by historical slip zones). These features together form a path distribution feature set, used as the master constraint for subsequent block response.

[0034] Based on the path distribution characteristics, the block deformation response is calculated and the mode type is determined. The block deformation response is based on the displacement field and settlement field: the global displacement field (obtained by finite element analysis and aggregated to the element scale, using a combination of element volume weighted average and key point displacement recording, with key points determined by the neighboring area of ​​the monitoring point, and the neighboring range determined by the upper limit of the monitoring point positioning error) is used to calculate the differential displacement between adjacent elements (obtained by decomposing the difference in displacement vectors on both sides of the shared boundary in the normal and tangential directions, with the decomposition direction determined by the interface normal direction). When the differential displacement exceeds the slip judgment threshold, it is marked as a significant relative slip area (the threshold is determined by the statistical upper limit of differential displacement under stable conditions and superimposed with a safety margin, which is determined by the upper limit of the monitoring system accuracy). The settlement field (extracted from the gravity component of the displacement field and distributed in a profile at key locations) is used to calculate the non-uniform settlement index (characterized by the settlement difference and settlement gradient between adjacent areas, with the settlement gradient obtained by the difference on the profile; the judgment threshold is taken as the abnormal boundary point between the settlement difference and gradient distribution, and verified by the allowable deformation control value of the key parts of the dam foundation). The cumulative effect along the main transmission channel is calculated (obtained by the cumulative displacement projection along the channel direction, which is determined by the main direction obtained from path tracing; the cumulative effect threshold is determined by the boundary point where the cumulative growth rate enters the abnormal interval in the construction stage sequence). The tensile stress control zone is identified (determined by the principal tensile stress being positive and continuously banded; the minimum continuous scale of the banded zone is determined by the structural characteristic scale; the tensile stress determination threshold is determined by the material tensile strength control value; the tensile strength control value is obtained from experimental values ​​after structural weakening reduction; the reduction range is determined by crack density and penetration index; the crack density and penetration index are obtained from crack grouping statistics). Under the constraints of the above response characteristics, the mode type is determined as follows: when the significant relative slip zone extends continuously along the same interface and highly coincides with the main control channel, it is determined to be a slip-dominant mode; when the tensile stress control zone extends continuously and coincides with the high-value area of ​​the concentration index, it is determined to be a tension cracking-dominant mode; when the compression-shear concentration zone forms a band-like extension along the channel and the uneven settlement index rises synchronously, it is determined to be a compression-shear coupling-dominant mode; when the upper displacement amplifies and the lower part forms a rotation center feature and the differential displacement shows a systematic sign change, it is determined to be a toppling or bending-dominant mode. The overlap index used in the above determination (determined by the spatial overlap ratio of the two types of areas, and the overlap ratio is obtained by grid or cell coverage counting) simultaneously determines the threshold when it appears. The threshold is taken as the low quantile control value of the overlap ratio of historical abnormal areas to ensure that the mode determination covers the real risk situation.

[0035] After determining the model type, the strain threshold corresponding to that model is obtained, and strain interaction analysis is performed on the rock mass unit to form a failure path sequence. The strain threshold is derived from dual constraints of materials and structure: the critical deformation control value of materials (determined by laboratory or in-situ tests, taking the unfavorable quantile value to reflect the influence of dispersion) and the structural weakening reduction (determined by a combination of fracture density, continuity, filling weakness, and water state; fracture density and continuity are obtained by fracture grouping statistics; filling weakness is determined by geological logging classification and controlled by the most unfavorable category; water state is determined by pore pressure and seepage data and taken under high water level control conditions) together form the strain threshold within the unit; the interface slip-related threshold is determined by the interface shear strength control value, which is jointly constrained by the interface frictional bond properties and the normal contact pressure level. The normal contact pressure level is obtained by the average pressure on the interface from the solution results, and the high quantile pressure is taken under unfavorable conditions to reflect the force transmission enhancement effect caused by unfavorable compaction. After thresholding is completed, strain interaction analysis is performed: First, at the mesh level, the strain distribution within the unit is extracted and the areas exceeding the strain threshold are marked. The marking process records the first time the limit is exceeded according to the construction stage or time sequence (the first time the limit is exceeded is determined by the stage number where the limit is first exceeded in the stage sequence); then, at the interface level, the relative strain state and relative slip trend of the two units are extracted to form the interface leading damage zone (the leading judgment threshold is the joint triggering condition of the interface slip correlation threshold and the strain threshold within the unit. The joint triggering requires that the two types of limits are spatially adjacent and temporally sequential. The temporal succession window length is determined by the loading... (Stage interval determined); then, the connectivity of the over-limit area is traced along the main control transmission channel. The minimum continuous scale for connectivity determination is determined (by the structural feature scale and converted to a uniform scale with the grid size, expressed as the number of grid layers, which is determined by dividing the scale by the grid size and taking the integer part), forming a sequence chain from the initial damage area to the through damage area; finally, the sequence chain is sorted based on the spatial advance distance in the channel direction at the first over-limit moment. The advance distance is obtained by accumulating the projection distance along the channel direction, thus forming a damage path sequence. The sequence simultaneously includes the spatial path and the evolutionary sequence.

[0036] Based on the failure path sequence, the stability indices of the block are integrated and the deformation patterns and failure paths of the geological blocks at the scale are output. The stability indices are integrated around the continuity index (determined by whether the failure path sequence connects the load-controlled area with the free surface or the exit area of ​​the weak zone, the exit area is determined by the superposition of boundary conditions and the geological weak zone), the control index (determined by the overlap ratio of the failure path sequence and the main control stress transmission channel, the overlap ratio threshold is taken as the low quantile control value of the overlap ratio under historical abnormal conditions), the expansion index (determined by the area increment and the advance distance increment of the over-limit area with the stage growth, the abnormal threshold of the increment is taken as the abnormal boundary point of the global increment distribution and checked by the acceleration change segment of the key part), and the consistency index (determined by the spatial overlap degree of the displacement concentration area, settlement concentration area and the monitoring abnormal area predicted by the model, the overlap degree threshold is jointly determined by the upper limit of the monitoring positioning error and the model grid resolution to ensure that the threshold matches the spatial resolution capability). When continuity, control, extensibility, and consistency all meet the threshold requirements, the final scale geological block deformation mode and failure path are output. The deformation mode is determined by the aforementioned mode determination results, and the failure path is the main sequence with the highest continuity and control among the failure path sequences. The location of its key control segment is marked simultaneously. The key control segment is determined by the overlap of the concentration index peak area, the interface leading damage area, and the monitoring anomaly core area. The threshold of the overlap segment length is determined by the structural characteristic scale to ensure that the key control segment corresponds to a continuous segment with engineering significance rather than a scattered point.

[0037] S3 includes obtaining threshold deviation detection results by comparing deformation modes with preset thresholds, performing iterative adjustment operations on the deviation data through grouping parameters to obtain an adjusted parameter set; for the adjusted parameter set, using a hierarchical clustering algorithm to restart the grouping process of rock mass unit data to generate preliminary rock mass unit groups; based on the preliminary rock mass unit groups, performing strain threshold acquisition operations from geological load data to determine the stress equilibrium adjustment scheme between units; using the stress equilibrium adjustment scheme between units, performing simulation input processing on the load data to determine changes in grid connection description; obtaining changes in grid connection description, integrating failure path sequences and stability indices through weighted calculations to obtain an optimized scale rock mass unit stiffness distribution map.

[0038] In this embodiment, the deformation mode is first compared with the preset threshold to obtain the threshold deviation detection result. The deformation mode is given by S2 and includes the main control mode type, the main control stress transmission channel, the displacement concentration area, the uneven settlement area, the interface slip risk area, and the failure path sequence. The preset threshold is used to limit the allowable response boundary of the project. Its components include the differential displacement threshold (the threshold is determined by the statistical upper limit of the differential displacement under stable working conditions, which is determined by the maximum fluctuation range of the stable time period, and superimposed with the upper limit of the monitoring system accuracy to form the final threshold, avoiding misjudging measurement noise as deviation), the uneven settlement threshold (the threshold is determined based on the allowable deformation control value of key parts, and checked with the statistical upper limit of the settlement difference under historical similar working conditions to ensure that the threshold is consistent with the project control requirements), and the stress concentration threshold (the threshold is taken as the abnormal boundary point of the distribution of the stress concentration degree index in the whole domain, and the abnormal boundary point is determined by the concentration degree index). The threshold is determined by the inflection point of the abrupt change from the concentrated interval to the tail, and is checked by the coverage rate of historical anomalies in key parts. If the coverage rate is insufficient, the threshold is adjusted to a more stringent direction. The connectivity threshold is determined by the proportion of the connectivity length of the disruptive path sequence to the length of the key influence area. The proportion threshold is taken as the low quantile control value of this proportion in historical instability or significant abnormal events to ensure that no risk is missed in the identification. The consistency threshold is determined by the spatial overlap ratio between the model-predicted anomaly area and the monitored anomaly area. The overlap ratio threshold is jointly constrained by the upper bound of the monitoring and positioning error and the model grid resolution to ensure that the threshold does not exceed the identifiable accuracy. The comparison process is carried out in priority order of key parts. First, key response quantities are extracted within the weak interlayer and fault influence zone of the dam foundation, and then key response quantities are extracted in the entire area. The exceedance magnitude of each key response quantity is calculated. The exceedance magnitude is characterized by the difference between the current response quantity and the corresponding threshold, and a deviation distribution is formed in space. The deviation distribution is clustered and identified. The clustering identification is constrained by spatial connectivity. The minimum connectivity scale (the scale is determined by the minimum scale of structural features, which is the minimum value of the thickness of the weak interlayer, the spacing between the main structural surfaces, and the width of the fault influence zone, and is converted into the number of grid layers as the criterion) is used to exclude isolated exceedance points. Finally, the threshold deviation detection results are formed. The detection results include the deviation location, deviation intensity, deviation coverage, and the overlap ratio between the deviation and the main control channel (the overlap ratio is obtained by counting the spatial overlap between the deviation area and the main control channel area. The count is based on the cumulative coverage area of ​​the unit or grid to avoid being dominated by individual grid points).

[0039] After obtaining the threshold deviation detection results, the deviation data is iteratively adjusted using grouping parameters to obtain the adjusted parameter set. Grouping parameters are used to control the hierarchical clustering grouping structure. They consist of feature weight parameters (the weights are determined by the feature set consisting of fracture orientation consistency, fracture density, continuity, fault influence intensity, and stiffness gradient intensity; the initial weights are obtained by ranking the correlation between the monitored anomalies in this area and the feature set; the correlation is determined by comparing the coverage of the anomaly area, and the higher the coverage, the higher the weight), spatial coherence constraint parameters (the parameters are determined by evaluating the spatial connectivity of the grouping results; the connectivity evaluation is measured by the proportion of the largest connected body formed by the same group of units; the proportion threshold is calculated by converting the smallest continuous unit scale required for engineering analysis), truncation level parameters (the parameters are determined by the transition point of the incremental difference within the group during the clustering and merging process; the transition point is determined by the inflection point where the increment enters the sudden increase interval from the stable interval, and is checked by the degree of boundary meshability; the check standard is that the boundary tortuosity does not exceed the mesh generation quality threshold, and the quality threshold is obtained by statistical analysis of the stable interval of mesh partitioning), and fault boundary constraint parameters (the parameters are determined by the contribution of the fault influence area to the deviation; the contribution is determined by the overlap ratio of the deviation area and the fault influence area; the higher the overlap ratio, the stronger the constraint). The iterative adjustment operation establishes directional correction rules driven by deviations: when the deviation is dominated by slip and the deviation area extends along the interface, the fault boundary constraint parameter is increased and the fault influence intensity characteristic weight is increased, so that the grouping results form clearer unit boundaries on both sides of the fault zone; when the deviation is dominated by stress concentration and the deviation area is located near the stiffness abrupt change zone, the stiffness gradient intensity characteristic weight is increased and the merging tendency across the abrupt change zone is reduced, so that the grouping results form a segmentation at the abrupt change zone; when the deviation is dominated by uneven settlement and the deviation area spans multiple unit boundaries, the spatial coherence constraint parameter is increased and the truncation level parameter is adjusted, so that the grouping scale matches the settlement control area scale. After each iteration, the deviation intensity decrease rate (obtained by the ratio of the difference between the total deviation intensity of the previous round and the total deviation intensity of the current round to the total deviation intensity of the previous round) and the deviation range contraction rate (obtained by the change ratio of the deviation coverage area) are recalculated. When the decrease rate and contraction rate enter the stable interval, convergence is determined and iteration stops. The threshold of the stable interval is the change amplitude corresponding to two consecutive rounds of change being lower than the upper limit of the monitoring noise, to avoid misjudging small changes caused by noise as effective convergence. If the decrease rate shows a reverse increase, it is rolled back to the parameter set of the previous round and the adjustment step size is reduced. The step size reduction ratio is determined by the reverse increase amplitude. The larger the increase amplitude, the more the step size is reduced, thereby avoiding over-correction that causes drastic fluctuations in the grouping structure.

[0040] For the adjusted parameter set, a hierarchical clustering algorithm is used to restart the grouping process of the rock mass unit data and generate preliminary rock mass unit groups. Restarting the grouping first involves reconstructing the similarity relationships of the rock mass unit data. The similarity is obtained by weighting the adjusted feature weight parameters. The weighting process uses normalized features to eliminate the influence of dimensions. The normalization range is determined by the upper and lower bounds of the feature statistics for the entire region. The upper and lower bounds are taken from the stable data range, and outlier extreme values ​​are removed. The threshold for removing outlier extreme values ​​is taken from the low and high quantiles of the feature distribution tail. Subsequently, a step-by-step merging process is performed, with the merging sequence based on the criterion of minimum increase in intra-group differences, until the truncation level parameters are met. After grouping, a geometric consistency check is performed: the number of spatial connected components is calculated for each group result. If the number of connected components exceeds a preset upper limit, fragmentation is determined. The fragmentation upper limit is calculated from the smallest continuous unit scale. Fragmented regions are merged. The merging principle is based on maximizing the boundary contact area with adjacent main connected components. The contact area is obtained by accumulating the shared boundary area, ensuring that the boundary is simpler and easier to mesh after merging. Minimal units are subjected to absorption processing. The threshold for determining minimal units is determined by the low quantile control value of the unit volume relative to the total volume of the region, and is checked against the minimum resolvable volume required for the mesh size to avoid generating units that cannot be stably partitioned.

[0041] Based on the preliminary rock mass unit grouping, strain thresholds are obtained from the geological load data, and an inter-unit stress equilibrium adjustment scheme is determined. Here, the strain thresholds are generated based on the load control conditions. Specifically, the geological load data is first decomposed into components: self-weight component, initial geostress component, water pressure component, construction stage component, and seismic inertia component. Then, the contribution of each component to key components is calculated. The contribution is obtained as the proportion of the increment of the key response caused by each component to the total increment. The increment is obtained by superimposing the components one by one on the same model and recording the response changes. The component combination with the highest contribution is selected as the control load combination. The strain threshold is determined under a controlled load combination. The threshold is obtained by reducing the critical deformation control value of the material. The critical deformation control value is determined by taking the unfavorable quantile value from indoor or in-situ tests. The quantile point of the unfavorable quantile value is determined by the data dispersion; the greater the dispersion, the more the quantile point is biased towards the unfavorable side. The reduction range is determined by the degree of structural weakening, which is obtained by comprehensively evaluating fracture density, continuity, and the degree of filling weakness. Fracture density and continuity are obtained by fracture grouping statistics. The degree of filling weakness is determined by geological logging classification and controlled by unfavorable categories. The water-bearing state is determined by pore pressure or seepage data and determined by high water level control conditions. After the strain threshold is determined, the strain imbalance index is calculated for the shared interface of adjacent units. The imbalance index is characterized by the difference in the degree of strain exceeding the threshold on both sides of the interface. When the difference exceeds the imbalance threshold, it is determined to be an imbalance interface. The imbalance threshold is determined by the statistical upper bound of the imbalance index under stable conditions and superimposed with the upper bound of the numerical solution error. The upper bound of the numerical solution error is determined by the upper bound of the difference between the results before and after mesh refinement. For the stress balance adjustment scheme between the units of the imbalance interface, the scheme includes the load components that need to be adjusted (determined by the component with the highest contribution), the adjustment area (determined by the range of the connected body of the imbalance interface, and the minimum scale for determining the connected body is obtained by converting the minimum scale of the structural features), and the adjustment amplitude level (the amplitude level is obtained by classifying the amplitude of the imbalance index exceeding the limit, and the classification threshold is taken as the upper boundary of the concentrated interval of the imbalance index distribution and the boundary point of the abnormal tail).

[0042] The load data is simulated and input processing is performed using an inter-element stress equalization adjustment scheme to determine changes in mesh connectivity. Simulation input processing is implemented in two aspects: spatial allocation and stage sequence. For areas requiring reduced local concentration, the load spatial allocation is reconstructed to create a smoother gradient between adjacent elements. The smoothness is constrained by the upper limit of the load gradient, which is determined by the anomalous boundary point of the gradient distribution under the control load combination. For areas requiring suppression of slip triggering, the load stage sequence is adjusted so that interface normal compaction is established before tangential drive enhancement. The magnitude of the stage sequence adjustment is constrained by the consistency between the construction stage timing and the inflection point of the monitoring response; if the consistency is insufficient, the adjustment is rolled back and the stage adjustment magnitude is reduced. After load reconstruction, the mesh connectivity description is resolved and compared: Geometric connectivity changes arise from changes in the relationships between adjacent elements due to regrouping. These changes are obtained by counting the differences between adjacent elements and the set. A significant change is determined when the difference count exceeds a threshold, controlled by the ratio of the number of elements in the previous grouping to the current grouping. The ratio control value is taken as the boundary point for entering the stable interval. Mechanical connectivity changes arise from changes in interface contact state and interface force transmission order. Changes in contact state are characterized by changes in the effective bearing area of ​​the interface, obtained by the ratio of contact elements to the total interface area. A change in ratio exceeding a threshold indicates a change in state, determined by the upper bound of proportional fluctuation under stable conditions plus the upper bound of numerical error. Changes in interface force transmission order are determined by the consistency of the ordering of the first few channels of the main control channel. The ordering consistency threshold is determined by the lower bound of the ordering stability statistics under historical stable conditions. If the threshold is lower than the lower bound, the main control channel structure is determined to have changed.

[0043] After confirming the changes in the mesh connectivity description, the failure path sequence and stability indices were integrated, and an optimized scale rock mass element stiffness distribution map was obtained through weighted calculation. The failure path sequence was obtained from the re-solved strain interaction analysis, and the sequence includes the initial damage segment, the extension segment, and the penetration segment, with each segment labeled with a stage number. The stability indices include penetration index, concentration index, slip risk index, settlement imbalance index, and consistency index. All index calculations were performed within the same mesh and element system to ensure comparability. The weights in the weighted calculation were determined simultaneously when the indices appeared: key location weight (the weight is determined by the coverage ratio of the key location, which is obtained by the area ratio of the high-value area of ​​the index within the key location; the higher the area ratio, the higher the weight), risk sensitivity weight (the weights of penetration and slip risk are higher than those of concentration and settlement imbalance; the weight ratio is determined by the proportion of the dominant mechanism in historical anomaly events; if the historical events are dominated by slip, the slip risk weight is increased), and consistency correction weight (when the consistency index is lower than the consistency threshold, the weight of the deviation-related indices is increased; the adjustment range is determined by the consistency deviation range; the larger the deviation, the greater the adjustment). The weighted results form a unit-level comprehensive evaluation value. Units with high comprehensive evaluation values ​​are considered regions that require improved resolution or adjustment of stiffness expression. When the comprehensive evaluation value is high and located in a section of the failure path, the stiffness zoning of the unit in that section is reconstructed. The zoning boundary is determined by superimposing the high stress concentration zone and the fault influence boundary, and is corrected by boundary meshability constraints to control the tortuosity of the zoning boundary. When the comprehensive evaluation value is high and the corresponding interface slip risk is prominent, the equivalent stiffness matching relationship of the units on both sides of the interface is re-evaluated. The re-evaluation is based on the interface force transmission strength and interface imbalance index, with the goal of making the strain exceeding the threshold on both sides of the interface tend to be balanced. When the comprehensive evaluation value is high and the corresponding settlement imbalance is prominent, the stiffness gradient of the settlement control area is reconstructed to make the gradient change consistent with the stratigraphic zoning and structural plane control direction. After completing the above unit-level adjustments, the stiffness distribution of the scale rock mass unit is regenerated, and a continuity check is performed across the entire region. The continuity check is based on the standard that the stiffness difference distribution of adjacent units enters the concentrated interval. The threshold of the concentrated interval is taken as the abnormal boundary point of the stiffness difference distribution in the whole region and checked with the coverage rate of the unfavorable area in the key part. After the check is passed, the optimized scale rock mass unit stiffness distribution map is output.

[0044] S4 includes obtaining material property parameters of the weak interlayer region of the dam foundation by optimizing the stiffness distribution map of the rock mass unit at the scale, thus obtaining local geological model data; constructing a mesh model using the finite element analysis method based on the local geological model data, and determining the boundary conditions and load distribution; inputting hydraulic seepage influence parameters from the mesh model, which are obtained in advance from the geological exploration data of the dam foundation, to calculate the stress field changes in the weak interlayer of the dam foundation. The calculation process uses the stress balance equation to integrate the seepage pressure distribution and obtain the local stress concentration value; based on the local stress concentration value, integrating the uneven settlement distribution data, if the settlement deviation exceeds the preset threshold, adjusting the model parameters to determine the settlement equilibrium state; obtaining the settlement equilibrium state, and combining it with the failure path simulation, determining the location of the potential sliding surface.

[0045] In this embodiment, the range of the weak interlayer region in the dam foundation is first determined based on the optimized scale rock mass unit stiffness distribution map, and material property parameters are obtained to form local geological model data. The range of the weak interlayer region is determined by the set of units corresponding to the interlayer geological identifiers in the optimized stiffness distribution map. The determination process is based on the spatial index of the unit boundary description, and the units belonging to the interlayer are aggregated for connectivity. Scattered and isolated units are eliminated. The minimum continuous scale for connectivity determination (which is obtained by converting the minimum value among the interlayer thickness, interlayer extension scale, and monitoring point spacing into the unit scale, and is used to exclude discrete fragments without engineering significance) is used to eliminate discrete fragments. After determining the extent of the interlayer, material property parameters are extracted. These parameters include at least stiffness parameters (obtained by pooling values ​​from the optimized stiffness distribution map within the interlayer units, using a volume-weighted approach where the unit volume is determined by the proportion of the unit volume to the total interlayer volume, thus maximizing the contribution of the main interlayer region to the representative value), density parameters (determined by indoor density tests corresponding to the interlayer lithology; if multiple batches of samples exist, weighting is applied based on the sample quantity and the coverage of the representative stratigraphic layer, with the coverage determined by the proportion of the sample's depth interval to the interlayer depth interval), and strength parameters (determined by indoor shear and triaxial tests matching the interlayer lithology; values ​​are preferentially selected from samples consistent with the interlayer's water-bearing state, with consistency confirmed by comparing water level conditions and pore water pressure monitoring records at the time of sampling; when significant test dispersion occurs, unfavorable quantile values ​​are used, with the quantile points determined by the degree of dispersion). The dispersion is measured by the ratio of the dispersion range of the test results to the median. The larger the ratio, the more the quantile point is biased towards the unfavorable side. Interface contact attribute parameters (determined by the roughness of the interface between the interlayer and the surrounding rock, the weakness of the filling material, and the aperture. The roughness is determined by the mapping and grading of the exposed surface. The weakness of the filling material is determined by the logging and classification. The aperture is determined by the statistics of borehole and cavern exposure and the representative value of the unfavorable side is taken. The representative value of the unfavorable side is determined by the high quantile control point of the aperture distribution to ensure that the interface shear resistance is not overestimated). Permeability-related parameters (pre-acquired and zoned by the geological exploration data of the dam site. The pre-acquired sources include pressure water test, permeability test, borehole water level observation and seepage interpretation results. The zone boundary is determined by the abrupt change position of the permeability difference in the test section. The abrupt change is judged by the change of permeability capacity of the adjacent test section entering the abnormal zone. The boundary point of the abnormal zone is determined by the tail inflection point of the statistical distribution of permeability capacity in the whole area. The above parameters together constitute the local geological model data. The local geological model data synchronously includes the material parameters of the surrounding rock units in the vicinity, with the interlayer as the center. The vicinity range is determined by the stress concentration influence range. The influence range is obtained by statistically analyzing the distance from the stress increment near the interlayer to the background level in the global model. The background level is determined by the stress mean of the same lithological region outside the interlayer.

[0046] After generating local geological model data, finite element analysis was used to construct a local mesh model and determine boundary conditions and load distribution. The local modeling scope extends outward from the interlayer geometric boundary to the boundary of the influence zone. The boundary of the influence zone is determined by boundary sensitivity verification constraints. The verification process involves establishing multiple sets of local models for different outward extension distances and comparing the magnitude of changes in key response quantities within the interlayer with the boundary position. Key response quantities include the maximum principal compressive stress, maximum shear stress, displacement gradient within the interlayer, and interface slip tendency intensity. When the change magnitude is lower than the change level corresponding to the upper limit of the monitoring resolution, the boundary position is deemed to meet the requirements. The upper limit of the monitoring resolution is jointly determined by the accuracy of the monitoring equipment and the upper limit of the fluctuation during the stable period, with a more stringent value being adopted. Mesh generation is constrained by the interlayer thickness direction and the interlayer boundary, ensuring that the interlayer thickness direction is discretized by at least multiple layers of elements. The number of multiple layers is obtained by rounding up the ratio of the interlayer thickness to the target mesh size. The target mesh size is determined by the smallest scale among the interlayer thickness, the pore pressure gradient change scale, and the stress concentration zone width, thereby ensuring that stress and pore pressure changes within the interlayer are captured. Local densification was implemented at the ends of interlayers, interlayer bends, and the interface between interlayers and high-stiffness surrounding rock. The extent of the densified zone was determined by superimposing the high-value stiffness gradient zone on the optimized stiffness distribution map with the continuous section of the failure path. The high-value stiffness gradient zone was determined by the position where the stiffness difference between adjacent elements entered the abnormal interval. The boundary point of the abnormal interval was determined by the tail inflection point of the statistical distribution of stiffness difference across the entire area. Mesh quality control was implemented throughout the meshing process. Quality indicators included minimum interior angle, element aspect ratio, and degree of geometric distortion. The threshold was determined by numerical solution stability verification. The verification method involved generating meshes for the same region using different quality thresholds and comparing the convergence stability and key response fluctuations. The most lenient threshold that entered the stable interval was selected as the final threshold to balance computational efficiency and stability. Boundary conditions and load distribution are determined in the local model according to the principle of global consistency: the far-field boundary displacement constraint is obtained by projecting the displacement field of the global model at the corresponding boundary. The projection process maps the nodal displacements on the global boundary to the local boundary nodes. The upper limit of the mapping error is controlled by the local mesh size. If the upper limit is exceeded, the local boundary is re-meshed to improve node alignment. The boundary equivalent stress load is obtained by extracting the stress field of the local outer boundary section of the global model and then performing area integration. The conversion process transforms the continuous distribution into local boundary nodal forces. The nodal force distribution is determined according to the boundary patch area ratio to conserve the boundary load. The self-weight load is calculated by the element density and gravity direction. The density is taken as the representative value of the element density in the aforementioned local geological model data. The dam body transmitted load is obtained by projecting the reaction force on the bottom surface of the dam body or the stress of the corresponding section of the global model. The projection area is determined by the dam foundation contact range. The contact range is determined by engineering geometry and contact surface exposure mapping. The water pressure-related boundary load is determined by the control water level. The control water level is determined by the high water level control point in the exploration water level statistics. The high water level control point is determined by the high quantile of the water level time series distribution to cover unfavorable working conditions.

[0047] After establishing the local mesh model and configuring the boundary loads, the hydraulic seepage influence parameters are introduced and the stress field changes within the interlayer are calculated to obtain the local stress concentration values. The parameters affecting hydraulic seepage include seepage zones, pore pressure boundaries, and seepage channel locations. Seepage zones are formed by segmented results of pressure water tests and seepage tests. The zone assignment uses segment representative values, which are the median values ​​of the test results within the same zone, and are purified using outlier removal rules. The outlier removal rules are constrained by the repeatability of the tests and the continuity of adjacent segments. Those with insufficient repeatability and differences from adjacent segments that fall into the outlier range are judged as outliers. The pore pressure boundary is jointly determined by the upstream control water level, the downstream water level, and drainage conditions. The upstream control water level is taken from the aforementioned high water level control point, and the downstream water level is taken from the statistically stable value during the operation period and superimposed with the upper limit of seasonal fluctuations. The upper limit of fluctuations is determined by the statistical upper limit of the stable period. The seepage channel locations are determined by superimposing the connectivity of geologically weak zones, fracture-dense zones, and weak interlayers. Fracture-dense zones are determined by the spatial distribution of high-density groups in the fracture grouping results, and the high-density discrimination threshold is determined by the tail inflection point of the fracture density statistical distribution. The calculation of the pore pressure field first involves solving for the seepage distribution in a local model based on the seepage zoning and water level boundaries, with the solution forming the pore pressure value within the unit. Subsequently, in the mechanical solution, the pore pressure is treated as a bulk force participating in the mechanical equilibrium, causing the stress state of the skeleton to adjust with changes in pore pressure. The mechanical solution outputs updated stress and displacement fields. If the permeability needs to be updated in conjunction with stress changes, the seepage zoning parameters are updated at the interlayer and its interface based on the concentrated shear deformation zone. The update magnitude is determined by the correlation evaluation between permeability and shear deformation. The correlation evaluation is determined by the degree of spatial overlap between historical seepage anomalies and deformation anomalies; the higher the degree of overlap, the larger the update magnitude. Then, the pore pressure field is re-solved and re-enters the mechanical solution, iterating until the changes in the pore pressure and displacement fields enter a stable range. The stable range criterion is jointly constrained by the change magnitude of multiple consecutive iterations being lower than the upper bound of the monitoring noise and the upper bound of the grid discretization error. The upper bound of the monitoring noise is obtained from the statistical analysis of monitoring fluctuations during the stable period, and the upper bound of the grid discretization error is obtained from the statistical analysis of the differences in key response quantities before and after grid refinement, with a more stringent value. Local stress concentration values ​​are extracted from the stress field after iterative convergence. The extraction process first searches for stress peak points inside the interlayer and near the interlayer interface. Then, a neighborhood is defined with the peak point as the center, and the representative stress level of the neighborhood is calculated. The neighborhood range is determined by the interlayer thickness and the mesh size. The representative stress level is taken as the median value of the neighborhood after removing the peak value to avoid the pull of single-point anomalies. The local stress concentration value is expressed as the ratio of the peak value to the representative value, and the mesh sensitivity is checked. The check method is to implement local densification in the peak region and then recalculate. If the peak position is stable and the ratio change enters the stable range, the concentration value is confirmed to be valid. The stable range is determined by the change before and after densification being lower than the upper bound of the mesh discretization error.

[0048] After obtaining local stress concentration values, the uneven settlement distribution data is integrated, and it is determined whether the settlement deviation exceeds a preset threshold. Then, a settlement equilibrium state is formed through parameter adjustments. The uneven settlement distribution data sources include monitored settlement data and settlement results calculated by local models. Before integration, the two are unified under a reference standard. The reference standardization converts the zero-point reference of different measuring points into the same reference surface. The reference surface is determined by the engineering control reference points. The control reference points are determined by the set of measuring points with the highest stability and far away from the interlayer influence zone. The stability judgment is based on the minimum long-term fluctuation amplitude. The calculation of settlement deviation revolves around key profiles, which are determined by the dam axis, the direction of the dam heel and toe, and the direction of the interlayer. The profile selection is based on the superposition result of the interlayer extension direction and the direction of the main control transmission channel. The direction of the main control transmission channel is determined by the stress transmission path extraction result of the previous step. On each profile, the settlement difference between adjacent measuring points, the location of settlement gradient abrupt change, and the degree of asymmetry of the settlement basin morphology are calculated. The settlement difference is obtained from the settlement difference between adjacent points. The settlement gradient abrupt change is identified by the abrupt increase area of ​​the settlement difference along the profile sequence. The abrupt increase discrimination threshold is determined by the statistical upper limit of the settlement difference during the stable period and superimposed with the upper limit of the monitoring accuracy. The degree of basin asymmetry is determined by the offset of the settlement center relative to the dam axis on the profile and the difference in settlement amplitude on both sides. The preset threshold is used to determine whether the settlement deviation has entered the risk zone. The determination of the preset threshold follows the dual constraints of engineering control and data statistics: first, the design-allowed uneven settlement control value is taken as the upper limit boundary, then the upper limit of the stable period monitoring fluctuation is taken as the lower limit boundary, and finally the threshold is taken as a value between the upper and lower limits that covers the minimum deviation level of historical abnormal events. The minimum deviation level of historical abnormal events is determined by the statistical low quantile of settlement deviation during historical abnormal periods, so that the threshold has risk coverage capability. When the settlement deviation exceeds the preset threshold, the model parameter adjustment is initiated and the settlement equilibrium state is determined. The parameter adjustment is carried out in order of sensitivity: First, the interlayer stiffness zoning parameters are adjusted (the zoning parameters are determined by the high value zone of the stiffness gradient inside the interlayer, the high value zone is determined by the position of the stiffness difference entering the abnormal interval, and the boundary point of the abnormal interval is determined by the inflection point of the tail of the stiffness difference statistics), and the adjustment direction is made so that the settlement gradient inside the interlayer is consistent with the geological zoning; then, the interlayer interface contact attribute parameters are adjusted (the contact attribute parameters are determined by the interface roughness classification, filling classification and aperture statistics), and the adjustment direction is made so that the intensity of the interface slippage trend decreases and is coordinated with the position of the stress concentration area; then, the pore pressure boundary and permeability zoning parameters are adjusted (the pore pressure boundary is determined by the control water level, and the permeability zoning is determined by the test segment), and the adjustment direction is made so that the high value area of ​​pore pressure is consistent with the position of the seepage channel revealed by exploration and suppresses the local abnormal peak value of pore pressure. After each round of parameter adjustment, the seepage and mechanical coupling solution is re-executed and the settlement deviation is recalculated. The determination of the settlement equilibrium state simultaneously meets two conditions: the settlement deviation falls back to within the preset threshold, and the change amplitude of the settlement deviation in multiple consecutive iterations enters the stable range. The stable range criterion is jointly determined by the upper bound of the monitoring noise and the upper bound of the grid discretization error, and a more stringent value is taken, so as to avoid misjudging numerical fluctuations as equilibrium.

[0049] After determining the settlement equilibrium state, the location of the potential sliding surface is determined by combining the failure path simulation. The failure path simulation results are derived from the failure path sequence output by the previous steps. The sequence is spatially mapped in the local model. The mapping process projects the through segment, extension segment and starting segment of the failure path onto the local mesh element set, and establishes a path neighborhood zone near the interlayer and its interface. The width of the neighborhood zone is jointly determined by the interlayer thickness and the width of the stress concentration zone. The location of potential sliding surfaces is determined by a combination of three conditions: The first condition is the force-driven condition, which is characterized by a continuous band-like region in space where high-value local stress concentrations are concentrated. The minimum continuous scale for determining the continuity of the high-value band is calculated by converting the interlayer extension scale and the grid size. The second condition is the seepage influence condition, which is characterized by the overlap of high-value pore pressure areas and high-shear deformation areas near interlayers or interfaces. The overlap is determined by the spatial overlap ratio, and the overlap ratio threshold is determined by the low quantile control point of the overlap ratio between historical seepage anomalies and deformation anomalies. The third condition is the evolution path condition, which is characterized by the failure path penetrating the above-mentioned high-value bands and forming a continuous chain from the loading control area to the free surface or weak zone outlet. The outlet area is jointly determined by the surface or cavern free surface and the interlayer exposure location. The connectivity of the set of elements satisfying the three conditions is traced and their geometric centerlines are extracted. The centerlines extend along the normal to form a planar band. The thickness of the planar band is determined by the interlayer thickness and the width of the interface influence. The width of the interface influence is obtained by the bandwidth statistics of the area with significant interface slip trend. The consistency of slip direction is further checked on the planar band. The consistency of direction is determined by the angle between the relative displacement direction of the element and the direction of the main control channel entering the concentration interval. The concentration interval is determined by the range of the main peak of the angle statistical distribution. After the check is completed, the position of the potential sliding surface is output. The position of the potential sliding surface is given in the form of a continuous spatial band, and the position of the control segment is marked simultaneously. The control segment is determined by the overlap of the stress concentration peak area, the pore pressure peak area and the failure path through section. The lower limit of the length of the overlap area is determined by the minimum scale of the structural features to ensure that the control segment has engineering significance of continuity.

[0050] S5 involves obtaining local stress concentration data through the deformation patterns of geological blocks at various scales, integrating relevant parameters of the deformation patterns, and determining the initial response distribution of the underground structure. For the initial response distribution, a multi-scale nested simulation algorithm is used to construct an overall mechanical response model by integrating relevant parameters of the deformation patterns and local stress data, resulting in a simplified macroscopic behavior distribution. Geological block-related information is extracted from the macroscopic behavior distribution to calculate groundwater level influence parameters and determine the preliminary location of the potential sliding surface. The preliminary location is obtained, and the relevant parameters of the deformation patterns are adjusted based on the local stress concentration data to determine the correction value for the sliding surface location. Based on the correction value, the relevant data of the macroscopic behavior distribution are integrated to obtain the prediction results of the macroscopic behavior of the underground structure.

[0051] In this embodiment, local stress concentration data is first obtained through the deformation model of a scale-scale geological block, and relevant parameters of the deformation model are integrated to determine the initial response distribution of the underground structure. The deformation model of the scale-scale geological block comes from the finite element solution and model classification output of the previous steps, and includes the dominant model type, dominant stress transmission channels, distribution of significant interfaces with relative slip between blocks, distribution of settlement imbalance, and distribution of failure path continuity segments. The local stress concentration data comes from the stress and pore pressure coupled solution output of the local model of key parts, including the stress peak zone in the weak interlayer, the force transmission abrupt change zone at the interface between the interlayer and the surrounding rock, the high pore pressure zone, and the shear deformation concentration zone. To achieve spatially consistent integration, the coordinate references and element indices of the global and local models are first aligned. This alignment is achieved through boundary description indices. Specifically, using the weak interlayer boundary and the fault influence zone boundary as anchor points, feature point sets on the boundaries are extracted from both models. These feature point sets consist of boundary inflection points and points with significant curvature changes. The selection threshold is taken as the abnormal boundary point of the boundary curvature statistical distribution, and the selection is checked with grid resolution to ensure that the number of selected points matches the model's resolution capability. Subsequently, the feature points are matched one by one, and the spatial mapping relationship is obtained so that the local stress concentration data can be projected onto the global element set and maintain positional consistency. After alignment, a set of parameters related to the deformation mode is formed. This set includes at least the following parameters: stiffness gradient parameters (derived from the stiffness difference between adjacent elements, where the stiffness difference is the difference in representative values ​​of the elements in the optimized stiffness distribution map; the anomaly detection threshold is the inflection point at the tail of the global stiffness difference distribution, used to mark abrupt change zones); connection strength parameters (derived from the force transmission strength at the interface between elements, which is obtained by integrating the interface reaction force along the normal and tangential directions, and normalized to obtain the unit area strength; normalization aims to eliminate the bias caused by differences in interface area); and interface slip tendency parameters (derived from the cumulative amount of relative displacement on both sides of the interface in the tangential direction, where the cumulative amount is...). The parameters are statistically analyzed according to the loading stage sequence. The threshold for significant trends is the upper limit of relative displacement fluctuation under stable conditions, superimposed with the upper limit of monitoring accuracy to avoid noise triggering. The parameters are: settlement gradient parameters (the parameters are obtained by the difference of the settlement field on the key profile, which is determined by superimposing the dam axis direction and the interlayer direction control direction; the threshold for sudden change is the upper limit of settlement gradient statistics during the stable period, superimposed with the upper limit of monitoring accuracy); and failure path stage parameters (the parameters are obtained by the stage number of the first exceedance of each segment in the failure path sequence, which is derived from the load stage division; the stage division boundary is determined by aligning the inflection point of the construction record and the monitoring response; and the threshold for inflection point discrimination is the upper limit of fluctuation in the stable section).The above parameters are superimposed with local stress concentration data to generate the initial response distribution of the underground structure according to spatial location. The generation rules for the initial response distribution are as follows: the high stress concentration zone, the zone with significant slip trend, the zone with abrupt settlement gradient change, and the section with continuous failure path are regarded as the core area of ​​the anomaly. The core area of ​​the anomaly is extended outward in the direction of stiffness gradient attenuation to form the anomaly influence zone. The extension distance is determined by the distance from the stress concentration value attenuation to the background level. The background level is taken as the average value of the concentration values ​​in the same lithological stable zone. The area outside the anomaly influence zone is classified as the normal response zone, thus obtaining the regionalized initial response distribution of the entire domain.

[0052] After obtaining the initial response distribution, a multi-scale nested simulation algorithm is employed to construct an overall mechanical response model and obtain a simplified macroscopic behavior distribution by integrating deformation mode-related parameters and local stress data. The multi-scale nested simulation algorithm consists of a global coarse-scale model and local fine-scale sub-models, which achieve closed-loop operation through boundary transfer and equivalent backpropagation. The material field of the global coarse-scale model is determined based on the optimized scale rock mass element stiffness distribution map. Boundary conditions and load paths are determined based on the overall working conditions. The load path includes self-weight, initial geostress, dam-transmitted stress, water pressure, and a phased working condition sequence. The boundaries of the phased sequence are determined by aligning monitoring inflection points with construction records. The selection of the local fine-scale sub-model is based on the initial response distribution, covering the key control sections of the anomaly core area and the anomaly influence area. The control sections are determined by the overlap of high-value stress concentration zones, high-value pore pressure zones, and the failure path continuity section. The overlap discrimination threshold is taken from the low quantile control point of the overlap ratio in historical anomaly events to ensure that no key control sections are missed. The global-to-local transmission adopts equivalent boundary constraints, which are calculated from the displacement and stress fields extracted from the global model on the local outer boundary section. The calculation process is distributed according to the area ratio of the boundary surface to ensure that the transmission amount is conserved and consistent with the boundary geometry. The local-to-global transmission adopts equivalent weakening parameter transmission. The equivalent weakening parameters include equivalent stiffness reduction factor (the factor is obtained by statistically analyzing the equivalent stiffness changes in the stress concentration and strain over-limit regions in the local model, and the statistics are aggregated with regional area weights to make the influence of the through section stronger), equivalent connection reduction factor (the factor is obtained by the interface force transmission attenuation ratio in the interface slip trend area, and the attenuation ratio is obtained by comparing the force transmission intensity per unit area of ​​the interface before and after slippage, and the time sequence is determined by the load stage number), and equivalent permeation influence factor (the factor is obtained by statistically analyzing the degree of reduction of effective force in the high pore pressure area, and the degree of reduction is obtained by the decrease of effective force in the high pore pressure area relative to the background area, and the background area takes the average value of the stable area). After applying the equivalent backpropagation factor to the corresponding region of the global model, the global response is resolved. This process of propagation and backpropagation is repeated until convergence. The convergence criterion is determined by the simultaneous entry of both global and local key response quantities into the stable interval. Global key response quantities include the intensity ranking of the main control channels and the variation amplitude of the maximum displacement and maximum settlement across the entire domain. Local key response quantities include the peak position and variation amplitude of local stress concentration values ​​and the variation amplitude of slip trend intensity. The stable interval threshold is determined by the combined constraints of the upper bound of the grid discretization error and the upper bound of the monitoring noise, with the more stringent value being selected. The upper bound of the grid discretization error is obtained by statistically analyzing the differences in results before and after grid refinement, and the upper bound of the monitoring noise is obtained by statistically analyzing the monitoring fluctuations during the stable period. After convergence, a simplified macroscopic behavior distribution is output. The macroscopic behavior distribution is based on the global displacement field, settlement field, distribution of the main control stress transmission channels, distribution of the main weakening zones, and distribution of slip candidate zones, while retaining the identification of key control segments.

[0053] After outputting the macroscopic behavior distribution, relevant information about geological blocks is obtained from the macroscopic behavior distribution, and groundwater level influence parameters are calculated to determine the preliminary location of potential sliding surfaces. The relevant information about geological blocks includes relative displacement zones between blocks (formed by spatially connecting the differential displacements of adjacent block units; the significant threshold for differential displacement is taken as the upper limit of differential displacement fluctuation under stable conditions, superimposed with the upper limit of monitoring accuracy), concentrated shear deformation zones between blocks (formed by connecting high-value areas of shear deformation index; the high-value threshold is taken as the inflection point at the tail of the shear deformation distribution and verified by the coverage rate of historical anomalies in key areas), normal compaction distribution between blocks (obtained from the statistical analysis of interface normal contact pressure, and normalized per unit area to form comparable indicators), and subsidence basin morphology (formed by extracting closed areas of subsidence field contour lines; the basin center is determined by the maximum subsidence point, and the basin boundary is determined by the boundary line where the subsidence contour lines enter the background area; the boundary threshold is taken as the average subsidence value of the stable area superimposed with the upper limit of stable fluctuation). The calculation of groundwater level influence parameters is organized around water level zoning, pore pressure distribution, and effective stress reduction zones: the control water level parameters are taken from the high water level control points of the water level time series during operation. The high water level control points are determined by the high quantile points of the water level distribution and verified with the highest water level traces revealed by exploration; the water level spatial zoning is determined by the upstream, downstream, and drainage zone boundaries. The boundaries are determined by the topography and drainage structure location and compared with the seepage interpretation results; the pore pressure distribution is jointly determined by the seepage zone and the water level boundary. The high pore pressure zone is identified by the inflection point threshold of the pore pressure distribution tail; the effective stress reduction zone is obtained by superimposing the high pore pressure zone and the stress field. The degree of reduction is determined by the decrease in effective stress relative to the background zone. The background zone is taken as the average value of the stable zone with the same lithology. By overlaying geological block information with groundwater level influence parameters, the preliminary location of potential sliding surfaces is identified. The initial judgment rules are as follows: select a strip-shaped area that meets the following criteria: concentrated and continuous shear deformation, significant reduction in effective stress, high overlap with the main control channel, and high continuity with the weakened zone as the preliminary location. The overlap ratio threshold is taken as the low quantile control point of the overlap ratio in historical anomaly events. The minimum continuity scale is obtained by converting the interlayer extension scale and the grid size, which is used to exclude scattered candidate segments. The preliminary location is then checked for directional consistency. The directional consistency is determined by the angle between the sliding direction of the candidate zone and the direction of the main control channel entering the concentrated interval. The concentrated interval is determined by the range of the main peak based on the angle.

[0054] After obtaining the initial location, the relevant parameters of the deformation mode are adjusted based on the local stress concentration data, and the correction value for the sliding surface position is determined. The correction aims at "positional offset, enhanced continuity, and convergence of control segments." First, the degree of overlap between the initial location and the local stress concentration zone is calculated. The degree of overlap is obtained by the ratio of the area of ​​their spatial overlap to the area of ​​the initial location. If the ratio is insufficient, correction is triggered. The trigger threshold is taken as the low quantile control point of this ratio in the historical anomalies of key parts. After triggering, the weakening weight in the connection strength parameter is increased to make the slip candidate zone converge towards the area of ​​significant interface force attenuation. The adjustment range of the weakening weight is determined by the size of the overlap gap; the larger the gap, the greater the adjustment. Then, the degree of overlap between the initial location and the high pore pressure zone is calculated. If the degree of overlap is high and the slip trend parameter is synchronously significant, the constraint weight of the equivalent permeability influence factor is increased to make the candidate zone form a more continuous connected chain in the pore pressure control zone. The increase in the continuity of the connected chain is characterized by the increase in the length of the connected body, which is obtained by the difference in the maximum connected length of the candidate zone before and after correction. Next, segments that cross high-stiffness zones at the initial location are eliminated and backfilled. The elimination criteria are that the stiffness gradient parameters of these segments enter an abnormally high value area and lack settlement gradient coordination. The abnormally high value threshold is taken as the inflection point at the tail of the stiffness gradient distribution. The coordination criterion is that the overlap ratio of the settlement gradient abrupt change zone and the shear deformation concentration zone reaches a threshold. The threshold is taken as the upper limit of the overlap ratio in the stable zone plus the risk margin. The above adjustments form the sliding surface position correction value. The correction value includes spatial offset (the offset is obtained by statistically analyzing the distance from the initial location centerline to the corrected centerline, with the median value taken and the maximum value on the unfavorable side marked), continuity enhancement (the enhancement is characterized by the increase in the maximum connected length and the decrease in the number of connected elements), and control segment position change (the change is determined by the spatial migration distance of the control segment overlap area, which is obtained by the difference between the center points of the control segments).

[0055] Finally, the macroscopic behavior distribution data are integrated based on the corrected values ​​to obtain the prediction results of the macroscopic behavior of underground structures. The integration process first writes the corrected slip surface position back into the slip candidate zone layer of the macroscopic behavior distribution to make the distribution of weakened zones at the macroscopic level consistent with the slip surface position; then, the sorting of the main control channel and weakened zone is updated according to the changes in the control segment position. The sorting update is determined by the channel strength and the degree of channel crossing the control segment. The degree of crossing is obtained by the ratio of the overlap length between the channel and the control segment to the channel length; then, the consistency of the displacement field and the settlement field is checked. The consistency check is based on the overlap ratio between the predicted anomaly area and the monitored anomaly area. The overlap ratio threshold is constrained by the upper bound of the monitoring positioning error and the macroscopic grid resolution. If the threshold is not reached, the adjustment is reversed and the weights related to consistency are increased. After the review is completed, the macroscopic behavior prediction results of the underground structure are output. The prediction results include the macroscopic distribution of displacement and settlement across the entire area, the distribution of the main stress transmission channels and the main weakening zones, the corrected location of the potential sliding surface and its control section, and the marking of the pore pressure control section, stress concentration control section and slip coordination control section in key parts, so that the macroscopic prediction results have the ability to provide engineering interpretation paths and risk location.

[0056] S6 includes obtaining relevant data on stiffness gradient distribution from geological structure prediction results; constructing a preliminary description of mechanical response law based on the stiffness gradient distribution and the changing characteristics of underground deformation trends, thus obtaining the initial distribution characteristics of the mechanical response law; obtaining the specific value range of slip coordination parameters based on the initial distribution characteristics of the mechanical response law; adjusting the boundary conditions of the mechanical response law by mapping the slip coordination parameters with local stress distribution, and determining the corrected response distribution range; integrating relevant information on macroscopic behavioral characteristics based on the corrected response distribution range; constructing an input dataset for the finite element analysis method under the influence of geological environmental variables, and obtaining preliminary comparison results of the finite element analysis method; obtaining the deviation value of the monitoring data comparison through the preliminary comparison results; adjusting the weight ratio of slip coordination parameters and stiffness gradient distribution according to the data matching accuracy requirements, and determining the distribution range of behavioral consistency values; if the distribution range of behavioral consistency values ​​meets the preset threshold, then combining the relevant data of local stress distribution and underground deformation trends to generate the final geological structure prediction adjustment scheme and determine the verification conclusion of macroscopic behavioral characteristics.

[0057] In this embodiment, stiffness gradient distribution data is first obtained from the geological structure prediction results. Then, a preliminary description of the mechanical response law is constructed based on the stiffness gradient distribution and the characteristics of underground deformation trends, yielding the initial distribution characteristics of the mechanical response law. The acquisition of stiffness gradient distribution data is based on the unit-level stiffness field, which is derived from the optimized stiffness distribution map and its spatial index in the macroscopic behavior prediction results. In specific implementation, firstly, all adjacent unit pairs in the domain are enumerated and their shared boundaries are extracted. The boundaries are uniquely determined by the unit boundary description. For each adjacent unit pair, the stiffness difference is calculated. The stiffness difference is obtained by subtracting the representative stiffness values ​​of the units on both sides. The representative stiffness values ​​are obtained by weighted aggregation of unit volumes. During aggregation, extreme sub-regions within the unit are removed. The threshold for extreme sub-regions is taken as the tail control point of the stiffness distribution within the unit and verified by the zone boundary revealed by exploration. Then, the boundary area of ​​the stiffness difference is normalized. The boundary area used for normalization is obtained by summing the areas of the shared boundary patches, so that abrupt changes with large contact areas contribute more to the gradient expression. Then, the normalized stiffness difference is spatially connected for tracing. The tracing is constrained by the connectivity of the shared boundary. The minimum connectivity scale is converted from the minimum structural feature scale in the entire domain. The structural feature scale is the minimum value of the thickness of the weak interlayer, the spacing between the main structural surfaces, and the width of the fault influence zone. In the conversion process, the scale is converted into the number of unit layers and rounded up to avoid discrete small segments being misjudged as control zones. The threshold for identifying high-value zones of stiffness gradient is determined simultaneously during the identification process. The threshold is taken as the inflection point where the statistical distribution of normalized stiffness difference across the entire domain transitions from the concentrated interval to the tail-end anomalous interval. This inflection point is determined by the location of a sudden increase in the slope of the distribution curve and is verified using the coverage rate of historical anomalous points in key areas. If the coverage rate is insufficient, the threshold is adjusted towards a stricter direction until the coverage rate reaches a preset lower limit, which is determined by the median coverage rate of historical anomalous points. Through the above calculations, the spatialized high-value zones of stiffness gradient distribution and their intensity levels are obtained.

[0058] The acquisition of underground deformation trend characteristics is based on the common scale expression of monitoring and prediction sequences. Specifically, the monitoring data is first synchronized in time and aligned to a benchmark. Time synchronization uses a unified sampling period and is aligned via interpolation. Benchmark alignment uses a stable set of benchmark points as a reference, determined by the measuring points with the smallest long-term fluctuation amplitude. The fluctuation amplitude threshold is taken from the low quantile control points of the stable period fluctuation statistics. Next, trend direction and trend rate are extracted from the monitored displacement and settlement sequences. The trend direction is obtained statistically from the main direction of the displacement vector in the plane, determined by the range of the main peak of the direction distribution. The trend rate is determined by the median value of the increments in adjacent time periods to weaken the influence of outliers. Then, trend change segments are identified, determined by the location of a continuous turning point in the trend rate or trend direction. The threshold for continuous turning points is the upper bound of the stable period noise superimposed with the upper bound of the monitoring accuracy, avoiding the identification of noise disturbances as change segments. These trend change segments are then projected onto the spatial profile and unit set to form trend change bands. After extracting the high-value zones and trend change zones of the stiffness gradient, a description of the mechanical response law is constructed. This description is organized using segmented spatial rules: when the trend change zone crosses the high-value zone of the stiffness gradient, the overlapping segment is identified as the response control segment, and the minimum length of the overlapping segment is calculated using the minimum structural characteristic scale; when the trend change zone coincides with a local stress concentration zone, the overlapping segment is identified as the force-driven segment, and the overlap threshold is taken as the low quantile control point of the overlap ratio in historical anomalous events; when the high-value zone of the stiffness gradient coincides with the boundary of a subsidence basin, the overlapping segment is identified as the settlement constraint segment. The boundary of the subsidence basin is determined by the boundary line where the settlement contour lines enter the background area, and the boundary threshold is taken as the average settlement value of the stable area superimposed with the upper limit of stable fluctuations. These identifications form the initial distribution characteristics of the mechanical response law, which at least include the spatial range and level of the control segment, the force-driven segment, and the settlement constraint segment.

[0059] Subsequently, based on the initial distribution characteristics of the mechanical response law, the specific value range of the slip compatibility parameters is obtained. The boundary conditions of the mechanical response law are adjusted through the correlation mapping between the slip compatibility parameters and the local stress distribution to determine the corrected response distribution range. The slip compatibility parameters describe the degree of deformation compatibility on both sides of the interface, and the parameter values ​​are obtained based on segmented statistical analysis of the interface. In specific implementation, the entire domain interface segments are enumerated according to the inter-element connection description, and for each interface segment, the differential displacement on both sides of the interface, the interface normal compression level, and the interface shear deformation level are extracted. The differential displacement is obtained by accumulating the difference in displacement between the nodes on both sides of the interface in the tangential direction. The accumulation is statistically calculated according to the load stage number, which is derived from the working condition stage division. The stage division boundary is determined by aligning the inflection points of the construction records and the monitoring response. The inflection point discrimination threshold is taken as the upper limit of the stable segment fluctuation. The normal compression level is determined by the average unit area of ​​the interface normal contact pressure, which is obtained by integrating the interface reaction area and dividing by the interface area. The shear deformation level is characterized by the relative tangential displacement gradient on both sides of the interface and the high-value zone of shear strain in the interface neighborhood. The neighborhood range is determined by the mesh size and the structural feature scale. Subsequently, the interface segments were classified according to their attributes. The attribute classification was determined based on the interface roughness level, the infill weakness level, and the water-bearing zone. The roughness and infill levels were determined by geological logging classification, and the water-bearing zone was determined by pore pressure distribution and control water level zoning. Within each attribute level, the distribution of slip coordination parameters was statistically analyzed, and the upper and lower bounds of the value range were determined. The upper bound was taken from the high quantile control point of the slip trend intensity distribution at that level. The quantile point of the high quantile control point was determined by the sample size and dispersion within the level; the greater the dispersion, the more the quantile point is biased towards the unfavorable side. The lower bound was taken from the upper control point of the slip trend intensity distribution at that level under stable conditions. The upper control point was determined by the maximum fluctuation range statistically analyzed during the stable period, superimposed with the upper bound of monitoring accuracy, ensuring that the lower bound could exclude normal fluctuations. Thus, the specific value range of the slip coordination parameters was obtained, with each value range corresponding one-to-one with the interface attribute level.

[0060] Local stress distribution is derived from the output of the local model of key parts, including stress concentration zones, interface force transmission abrupt change zones, and shear deformation concentration zones. When establishing the correlation mapping between slip compatibility parameters and local stress distribution, the high-value zones of the local model are first projected onto the global interface segment set. This projection is achieved through boundary anchor point alignment. The anchor points are selected from the feature point set of the weak interlayer boundary and the fault influence zone boundary. The selection threshold for the feature point set is the boundary curvature anomaly boundary point, verified by the mesh resolution. Next, for each interface segment, the slip compatibility parameters are calculated to see if they fall into the unfavorable value range, and simultaneously, it is determined whether the interface segment is within the coverage range of the local high-value stress zone. When the overlap ratio between unfavorable slip segments and high-value stress segments exceeds the mapping threshold, a strong correlation is established. The mapping threshold is taken from the low quantile control point of the overlap ratio in historical anomaly events, verified by the mesh resolution to ensure that the threshold is not lower than the minimum identifiable overlap area ratio. Based on strong correlation, the boundary conditions of the mechanical response law are modified. The modification includes boundary constraint allocation and interface constraint attribute adjustment: when the strongly correlated segment is concentrated in the boundary neighborhood of the response control segment, the boundary constraint allocation is adjusted to reduce the constraint stiffness gradient in that neighborhood, so that the boundary force transmission path is consistent with the control segment. The adjustment range is determined by the coverage ratio of the strongly correlated segment; the higher the coverage ratio, the more adjustment is made. When the strongly correlated segment runs through the interface, the interface constraint attribute is adjusted to ensure that the tangential force transmission capacity remains consistent with the change in normal compression level. The adjustment direction is determined by the relative position of the slip coordination parameter within the value range; the closer it is to the upper bound, the more emphasis is placed on the weakened expression. After the modification is completed, the modified response distribution range is formed, which expands outward from the strongly correlated control segment as the core. The expansion distance is the distance from the attenuation of the stress concentration value to the background level. The background level is the median of the concentration values ​​in the stable zone of the same lithology. The attenuation discrimination threshold is the position where the concentration value exceeds the background level and enters the stable noise interval. The noise interval is jointly determined by the upper bound of the grid discretization error and the upper bound of the monitoring noise.

[0061] After obtaining the corrected response distribution range, relevant information on macroscopic behavioral characteristics was integrated and a finite element input dataset was constructed under the influence of geological environmental variables to obtain preliminary finite element comparison results. The relevant information on macroscopic behavioral characteristics was extracted from the macroscopic prediction results and included displacement field, settlement field, main control channel, weakened zone, potential slip surface candidate zone, and control segment identifier. Geological environmental variables included at least groundwater level zoning, pore pressure distribution, permeability zoning, temperature zoning, and load paths during the construction phase. Groundwater level zoning is obtained from the time-series statistics of groundwater levels during operation. The upper limit is taken from the high water level control point, which is determined by the high quantile of the water level distribution and verified with the highest water level traces revealed by exploration. The lower limit is taken from the stable control point, which is determined by the midpoint of the water level distribution and compared with the normal water level during operation. Pore pressure distribution is determined by the water level boundary and the permeability zone. The permeability zone is formed by the segmented statistics of water pressure test and permeability test. The segment boundary is determined by the location where the permeability suddenly enters the abnormal zone. The boundary point of the abnormal zone is taken as the inflection point at the tail of the permeability statistics. Temperature zoning is obtained from the time-series statistics of temperature monitoring points. The unfavorable temperature difference control value is taken and the zoning is formed by spatial interpolation. The interpolation range is determined by the coverage of the monitoring points. Areas with insufficient coverage are supplemented by the representative values ​​of the same lithology zoning and marked with low confidence. When constructing the finite element input dataset, the computational domain is first trimmed with the corrected response distribution range, concentrating computational resources on the control section and its influence zone. Then, the material partitions are refined based on the high-value bands of the stiffness gradient, ensuring independent assignment on both sides of the abrupt change zone. Next, the interface attribute partitions are refined based on strongly correlated interface segments, giving the slip-incompatible segments independent constraint attributes. Boundary constraints are based on the corrected boundary condition allocation. Load combinations are taken from the working condition combinations corresponding to the monitoring period, with the working condition combination sequence aligned with the construction stage. The upper limit of the stage alignment error is taken as the upper bound of the monitoring sampling period, and verified by the inflection point identification accuracy. The finite element solution outputs preliminary comparison results, including predicted displacement, predicted settlement, interface slip trend, and the location and intensity of stress concentration zones. Resampling is performed according to monitoring points and monitoring profiles to ensure a comparison with the monitoring data of the same caliber.

[0062] Based on the preliminary comparison results, the deviation values ​​of the monitoring data were obtained, and the weight ratios of the slip coordination parameters and stiffness gradient distribution were adjusted according to the matching accuracy requirements to determine the distribution range of the behavior consistency values. The deviation values ​​were calculated item by item according to the indicators: the displacement deviation was calculated as the absolute difference at the monitoring point and normalized to form the relative deviation according to the representative value of the monitoring magnitude. The representative value was the median response of the monitoring point in the corresponding working condition period; the settlement deviation was calculated as the difference of the settlement difference curve on the key profile, and the maximum difference and the location of the difference were extracted. The location deviation was measured by the profile distance; the interface slip deviation was divided into directional deviation and strength deviation. The directional deviation was obtained by statistically analyzing the angle between the slip directions. The angle concentration interval was determined by the range of the main peak of the angle distribution. The strength deviation was obtained by the difference in slip trend strength and normalized according to the upper limit of the fluctuation during the stable period; the stress concentration location deviation was obtained by statistically analyzing the distance between the center lines of the high-value zone. The median value of the distance statistics was taken and the maximum value on the unfavorable side was marked. The matching accuracy requirement is jointly determined by the monitoring system accuracy, the allowable deviation of the project, and the risk identification requirements. The monitoring system accuracy is given by the equipment calibration and verified by the noise during the stable period. The allowable deviation of the project is determined by the design control value. The risk identification requirements impose stricter constraints on the key control sections. The degree of strictness is determined by the sensitivity of historical anomalies to the deviation. The sensitivity is obtained by statistically analyzing the slope of the anomaly identification coverage as a function of the deviation threshold. The larger the slope, the stricter the matching requirements of the key section. The weight ratio adjustment is driven by the sensitivity of the deviation source: when the deviation is mainly concentrated in the interface slip direction and intensity, the weight of the slip coordination parameter is increased and the weight of the stiffness gradient is decreased. The adjustment range is determined by the grading of the slip deviation exceeding the limit. The grading threshold is taken as the boundary point where the slip deviation distribution enters the tail of the anomaly. When the deviation is mainly concentrated in the stress concentration location and the morphology of the settlement basin, the weight of the stiffness gradient is increased and the slip coordination weight is suppressed. The suppression range is determined by the grading of the deviation exceeding the limit in the settlement and concentration location. The grading threshold is also taken as the inflection point of the tail of the deviation distribution. The behavioral consistency value is obtained by weighted aggregation of multiple deviations. The aggregation process first normalizes all deviations to the same scale, with the upper limit of the normalization scale being the allowable deviation for engineering. Beyond the upper limit, the deviation is truncated according to the unfavorable side to avoid extreme outliers dominating consistency. Then, the consistency value is formed by summing according to the weights. The distribution range of the consistency value is determined by the low quantile control point and the high quantile control point of the overall consistency value. The low quantile control point is used to express the performance of the unfavorable side, and the high quantile control point is used to express the overall level. The quantile points of the two are determined by the weight ratio of key parts. The higher the weight of the key parts, the more the low quantile points are biased towards the unfavorable side.

[0063] When the distribution range of behavioral consistency values ​​meets the preset threshold, the final geological structure prediction and adjustment scheme is generated by combining local stress distribution and underground deformation trend data, and the macroscopic behavioral characteristic verification conclusion is determined. The preset threshold is used to determine whether consistency meets the standard, and its determination follows the multi-source error synthesis constraint: the lower limit is jointly determined by the upper limit of monitoring noise and the upper limit of grid discrete error, and the more stringent value is taken. The upper limit of monitoring noise is obtained from the statistical analysis of monitoring fluctuations during the stable period, and the upper limit of grid discrete error is obtained from the statistical analysis of the difference in key responses before and after grid densification; the upper limit is determined by the engineering control accuracy requirements, and the key control section adopts a tightened threshold. The tightening magnitude is determined by the slope of the sensitivity of historical abnormal events, and the larger the slope, the stronger the tightening. When determining compliance, it is simultaneously checked that the consistency value of the key control section is not lower than the tightened threshold, and it is also checked that the low quantile control points in the distribution range of the overall consistency value meet the threshold requirements to avoid overall compliance but mismatch in key sections. After achieving the target, a final adjustment plan is generated. This plan is output in an itemized structure, with each item including at least the following: stiffness gradient zoning correction item (the correction targets are the unit zoning boundaries and representative stiffness near high stiffness gradient zones, based on stress concentration location deviation and subsidence basin deviation); slip compatibility parameter value range correction item (the correction targets are the upper and lower bounds and grading rules of strongly correlated interface segments, based on slip direction deviation and strength deviation); and boundary condition allocation correction item (the correction targets are the constraint allocation and equivalent load allocation of the outer boundary of the response distribution range, based on the consistency of the main control channel traversing the control section). Simultaneously, a macroscopic behavioral characteristic verification conclusion is generated. This conclusion clearly indicates that the predicted results and monitoring data meet the consistency requirements within the key control section and the entire domain, and clarifies that the synergistic relationship between potential slip surface candidate zones, stress concentration zones, high pore pressure zones, and trend change zones satisfies the control logic.

[0064] S7 includes obtaining specific data on fault reconstruction correlation from fracture grouping iterations. Fracture grouping iterations involve multiple iterative classifications using the geometric distribution and stress field data of scale-based fractures. Fault reconstruction correlations are calculated using interlayer displacement vectors and friction coefficients to determine correlation strength. Based on the results of matching degree threshold verification, an initial framework for feedback loop optimization is constructed to obtain a preliminary distribution for model consistency adjustment. For this preliminary distribution, relevant information from scale feature simulations is integrated. Scale feature simulations use the finite difference method to process fracture density and permeability data. An interlayer interaction correction process is then implemented, adjusting the interaction vector based on displacement gradient and shear modulus. The boundary conditions for the deviation correction process are used to determine the correction range for integrated geological evaluation. Based on the correction range, the input dataset for regional model construction is obtained, which includes the corrected stress distribution and deformation modulus. Under the influence of the coordination parameter fusion, which combines stiffness gradient and slip parameters through a weighted average method, deformation trend calibration is performed to obtain the final scheme optimization comparison results. By comparing the results, the range of behavioral consistency values ​​is determined. If it is lower than the preset threshold, the fracture grouping iteration and fault reconstruction association are re-executed, which includes updating the stress field data and friction coefficient to obtain the final regional geological evaluation model.

[0065] In this embodiment, specific data related to fault reconstruction are first obtained from the fracture grouping iteration. The fracture grouping iteration uses scale fracture geometric distribution and stress field data as joint inputs. The scale fracture geometric distribution data includes fracture spatial location, strike and dip characteristics, dip angle characteristics, spacing density, extension scale, and connectivity. This geometric distribution data is formed by field logging, borehole interpretation, and exposed surface mapping, and is then transformed into the same spatial framework after unified coordinate reference. The stress field data is obtained by solving the preceding macroscopic model or local model, and at least includes the main control force transmission direction, interface normal compaction level, interface tangential driving level, stress concentration zone location, and intensity ranking. Before iteration, scale unification is performed on the geometric and stress features. Scale unification is achieved by establishing a global statistical range for each feature. The upper and lower bounds of the statistical range are taken from stable data intervals, and extreme outliers are removed. The outlier removal threshold is taken from the inflection point at the tail of the feature distribution and checked against the upper limit of the acquisition error. The upper limit of the acquisition error is jointly determined by the mapping accuracy and the borehole positioning error, and a more stringent value is taken. After scale unification, a comprehensive similarity is calculated for each fracture or fracture set. The comprehensive similarity is constrained by both geometric similarity and stress similarity. Geometric similarity is formed by azimuth difference, spacing difference, extension scale difference, and spatial proximity. The azimuth difference is calculated using a directional continuity rule to ensure that the difference across directional boundaries remains continuous. Spatial proximity is determined by both three-dimensional distance and stratigraphic constraints. Stratigraphic constraints are determined by stratigraphic interface index to prevent fractures from different stratigraphic levels from being mistakenly merged due to similar projections. Stress similarity is formed by the stress level, shear drive intensity, and normal compaction level of the fracture neighborhood. The neighborhood range is determined by both the grid size and the minimum scale of structural features. The minimum scale of structural features is the minimum value of the weak interlayer thickness, the spacing between major structural planes, and the width of the fault influence zone, which is then converted into the number of grid layers. Subsequently, multiple rounds of iterative classification are performed. Each round of classification forms fracture groups based on maximizing comprehensive similarity. After each round, the stability of the groups is calculated. Group stability is jointly measured by the change ratio of group members and the change ratio of group boundaries. The change ratio of members is obtained by the crossover ratio of fractures in the same group in the current round and the previous round, and the change ratio of boundaries is obtained by the area difference ratio of the spatial envelope boundaries of the same group. The stability criterion threshold is determined simultaneously in the criterion. The threshold is the allowable change level after combining the upper bound of fracture acquisition error and the upper bound of grid discretization error. The upper bound of grid discretization error is obtained by statistically analyzing the differences in key stress and displacement indicators before and after grid refinement, and the composite value is the most stringent value. If the stability criterion is met for two consecutive rounds, the iteration terminates and the fracture grouping iteration results are output. At the same time, the proportion of spatial connectivity in each group is output. The proportion of connectivity is obtained by the ratio of the number of fractures with the largest connectivity to the total number of fractures in the group, which is used to determine the reliability of subsequent fault reconstruction and seepage channel identification.

[0066] During the fracture grouping iteration process, the fault reconstruction correlation is calculated synchronously to form correlation strength data. The fault reconstruction correlation uses interlayer displacement vector and friction coefficient as core inputs. The interlayer displacement vector is formed by the displacement difference between the upper and lower layers or the blocks on both sides at the same spatial position. The displacement difference is taken as the difference between the monitored displacement sequence or the numerical displacement field at the same position. Time synchronization and noise suppression are performed before calculation. Time synchronization is aligned by a unified sampling period. Noise suppression uses the upper limit of the stable period fluctuation as the filtering scale constraint. The upper limit of the stable period fluctuation is obtained by statistically analyzing the maximum fluctuation range of the displacement difference during the stable period and superimposing the upper limit of the monitoring accuracy. The friction coefficient is determined based on the interface roughness level, the filling weakness level, and the water-bearing state zoning. The roughness level and filling level are obtained by geological logging classification. The water-bearing state zoning is obtained by pore pressure distribution and control water level zoning. The friction coefficient value adopts the representative value of the unfavorable side. The representative value of the unfavorable side is determined by the low quantile control point of the friction value distribution of the same interface. The quantile of the low quantile control point is determined by the sample size and the degree of dispersion. The greater the dispersion, the more the quantile point is biased towards the unfavorable side. The calculation of correlation strength is carried out in segments: First, the directional consistency of inter-layer displacement vectors is judged. Directional consistency is determined by the angle between adjacent displacement vectors entering the concentration interval, which is determined by the range of the main peak of the angle distribution. Next, the synergy of displacement amplitude is judged. Synergy is determined by the correlation degree of displacement amplitude increments within the strip region. The correlation degree threshold is the upper bound of the correlation degree during the stable period plus a noise margin, which is determined by the upper bound of the time synchronization error. Finally, the friction coefficient is introduced to correct the slip tendency. The correction rule is to increase the correlation strength level in the low friction zone and decrease the correlation strength level in the high friction zone. The adjustment range is determined by the quantile position of the friction coefficient in its global distribution. The more unfavorable the quantile, the stronger the adjustment. Through the above calculations, the distribution of fault reconstruction correlation strength, the location of strong correlation zones, and the continuity index of strong correlation zones are obtained. The continuity index is characterized by the maximum connectivity length and the number of connected elements in the strong correlation zone. The maximum connectivity length is obtained by the projection length of the strip region along the main direction, which is determined by the main peak of the displacement vector direction.

[0067] Based on the matching degree threshold verification results, an initial feedback loop optimization framework was constructed, and a preliminary distribution for model consistency adjustment was obtained. The matching degree threshold is used to determine the consistency between the model prediction and the monitoring response. Its determination adopts a triple constraint: the upper bound of monitoring noise (obtained from the statistical analysis of monitoring fluctuations during the stable period and superimposed with the upper bound of monitoring accuracy), the upper bound of numerical discrete error (obtained from the statistical analysis of the differences in key indicators before and after grid densification), and the engineering control accuracy requirements (determined by the design control indicators). The strictest of the three is taken as the threshold benchmark, and a tightening coefficient is applied to key parts. The tightening coefficient is determined by the sensitivity of historical abnormal events to errors. The sensitivity is obtained by statistically analyzing the slope of the anomaly identification coverage rate as a function of the threshold. The larger the slope, the stronger the tightening. The matching degree verification was implemented using a consistency index system. The consistency indexes included the displacement field overlap ratio, subsidence basin location offset, slip trend direction deviation, and stress concentration zone location deviation. The displacement field overlap ratio was obtained by the ratio of the overlap area of ​​the predicted anomaly area to the union area of ​​the monitored anomaly area; the anomaly area threshold was formed by superimposing the background mean of the stable area with the upper limit of fluctuations. The basin location offset was obtained by the distance from the subsidence center point, with the center point being the point of maximum subsidence. The slip trend direction deviation was obtained by the median value of the slip direction angle and the maximum value on the unfavorable side. The stress concentration zone location deviation was obtained by the median value of the distance from the centerline of the high-value zone and the maximum value on the unfavorable side. If the matching degree was lower than the threshold, a feedback loop optimization initial framework was constructed. The framework superimposed the fracture grouping change sensitive area (the sensitive area was formed by the concentrated area of ​​changes in grouping members in two consecutive rounds, with the concentrated area threshold being the high quantile control point of the stability criterion threshold), the fault strong correlation zone, the stress concentration zone, the high seepage influence zone, and the monitoring deviation accumulation zone to form a preliminary distribution for consistency adjustment. The high-impact zone of seepage is formed by the overlap of high-permeability and high-pore pressure zones, with the high-value threshold taken as the inflection point at the tail end of the permeability and pore pressure distribution. The monitoring deviation cluster zone is formed by the spatial connectivity of deviation exceeding limits, with the minimum scale of connectivity criterion being a combination of the monitoring point spacing and the minimum scale of structural characteristics, with a more stringent value selected. The initial distribution is sorted by deviation intensity, which is obtained by weighting the exceedance magnitude and the coverage area. The coverage area is calculated by accumulating the unit area, giving priority to deviation zones with larger areas in the sorting.

[0068] After obtaining the preliminary distribution, relevant information from scale feature simulation was integrated, and the integrated correction range for geological evaluation was determined. The scale feature simulation employed the finite difference method to process fracture density and permeability data. First, a differential grid aligned with the preliminary distribution was established. The grid scale was constrained by the bandwidth of the preliminary distribution and the width of the seepage channels, with a more stringent scale chosen to ensure the differentiation of zonal features. Subsequently, fracture density and permeability were projected onto the differential grid. The projection used area-weighted aggregation, with the area weight determined by the overlap ratio between the grid cell and the area covered by the original data. Anomaly cleansing was performed on the projection results. The anomaly cleansing threshold was taken as the inflection point at the tail of the statistical distribution and verified against the location of high-permeability channels revealed by exploration. If no coverage was found, the threshold was adjusted towards a more stringent direction. After projection, fracture density and permeability gradients were calculated on the differential grid. The gradient direction was used to identify channel orientation, and gradient abrupt changes were used to identify boundary segments. The abrupt change threshold was taken as the upper bound of the gradient distribution in the stable zone, superimposed with the upper bound of the discrete error. The process then proceeds to interlayer interaction correction, which adjusts the interaction vector based on displacement gradient and shear modulus. The displacement gradient is obtained by differentiating the interlayer displacement differences on a differential grid and smoothed within a noise upper bound, which is statistically derived from the fluctuations in displacement differences during the stable period. The shear modulus is determined by consistency screening through material testing and inversion, prioritizing samples with the same lithology, water content, and stress level. If these conditions are not met, the representative value of the unfavorable side is used, determined by low-quantile control points in the distribution of similar samples. The direction of the interaction vector is taken from the principal direction of the displacement gradient, and the intensity of the interaction vector is expressed as the grade of the displacement gradient intensity converted to shear modulus. The grade classification threshold is taken from the inflection point at the tail of the interaction intensity distribution. The boundary conditions for the deviation correction process are adjusted based on the interaction vector. This adjustment is implemented at the outer boundary of the initial distribution, ensuring that the force transmission direction at the boundary aligns with the principal direction of the interaction vector. The adjustment range is determined by the deviation concentration intensity, which is composed of the deviation coverage area and the exceedance magnitude. After boundary adjustment, a corrected value range is formed. The corrected value range is used to limit the range of subsequent regional model parameter updates. The upper and lower bounds are jointly constrained by the historical stable operating condition statistical interval and the unfavorable operating condition statistical interval. The stable interval is determined by the parameter fluctuation range during the stable period, and the unfavorable interval is determined by the parameter inversion range in historical abnormal events. Finally, the overlapping area of ​​the two is taken as the usable corrected value range to avoid over-extrapolation.

[0069] The input dataset is constructed based on the corrected value range of the regional model, and the final optimization comparison results are obtained by performing coordinated parameter fusion and deformation trend calibration. The input dataset contains the corrected stress distribution and deformation modulus. The stress distribution is obtained by solving the updated global mechanics problem. The updated solution uses the results of strong correlation bands and interaction vectors as constraints to ensure that the path of the main control channel through the strong correlation band remains continuous. The deformation modulus is updated within the corrected value range, and the update magnitude is determined by the deviation intensity ranking. The higher the ranking, the larger the update magnitude. After the update, a continuity check is performed. The continuity check is based on the standard that the modulus difference distribution of adjacent elements enters the concentrated interval. The threshold of the concentrated interval is the inflection point when the modulus difference distribution enters the tail abnormal interval from the concentrated interval. The coordination parameter fusion is performed synchronously when constructing the input set. The fusion adopts a weighted average method to combine the stiffness gradient and slip parameters: the stiffness gradient is the value of the stiffness difference between adjacent elements after normalization by the boundary area; the slip parameter is the comprehensive level of the interface slip trend intensity and the synergy of the inter-layer displacement vector, and the synergy is determined by both directional consistency and amplitude synergy; the weight of the weighted average is allocated according to the source of deviation. If the deviation is mainly reflected in the positional shift of the subsidence basin and the positional shift of the stress concentration zone, the weight of the stiffness gradient is increased. If the deviation is mainly reflected in the mismatch between slip direction and slip intensity, the weight of the slip parameter is increased. The weight adjustment range is determined by the deviation exceeding the limit classification. The classification threshold is taken as the inflection point of the tail of the deviation distribution and is tightened for key parts. After completing the coordination parameter fusion, deformation trend calibration is performed. The calibration process targets the monitored trend change zone and compares the overlap ratio between the predicted and monitored trend change zones. The overlap ratio is obtained by the ratio of the overlap area to the union area. If the overlap ratio is insufficient, the modulus zoning and interface weakening level are adjusted within the correction range. The adjustment direction is to make the trend change zone converge towards the strongly correlated zone and the high-influence area of ​​seepage. The convergence criterion is that the overlap ratio increases for two consecutive rounds and enters the stable range. The stable range threshold is jointly determined by the upper bound of the monitoring noise and the upper bound of the numerical discretization error. After calibration, the final scheme optimization comparison results are output. The comparison results include the updated displacement, settlement, slip trend, stress concentration zone, and trend change zone, and are output in the form of monitoring points aligned with key profiles.

[0070] The behavioral consistency value range is determined by comparing the results and triggering cyclic updates or outputting the final model. The behavioral consistency value is obtained by weighted aggregation of multiple deviations, including displacement deviation, settlement deviation, slip direction deviation, slip strength deviation, and stress concentration location deviation. The upper limit of the deviation normalization is taken as the engineering allowable deviation, and deviations outside the upper limit are truncated according to the unfavorable side to prevent extreme values ​​from dominating. The weights are integrated with the aforementioned coordination parameters to ensure that the deviation aggregation and parameter update directions are consistent. The behavioral consistency value range is determined by the low quantile control point and the high quantile control point of the global consistency value. The low quantile control point is used to express the level of the unfavorable side, and the high quantile control point is used to express the overall level. The position of the quantile point is determined by the weight ratio of key parts. The preset threshold is used to determine whether it is necessary to re-execute the crack grouping iteration and fault reconstruction correlation. The preset threshold is determined by the strict side principle, taking the most stringent one among the upper limit of monitoring noise, the upper limit of numerical discrete error, and the engineering control accuracy requirement as the threshold benchmark, and applying a tightened threshold to the strongly correlated control segment. The tightening magnitude is determined by the slope of historical anomaly sensitivity. If the behavior consistency value range is lower than the preset threshold, the fracture grouping iteration and fault reconstruction correlation are re-executed. During the re-execution, the stress field data and friction coefficient are updated. The stress field data is obtained by resolving the latest modulus zoning, interface weakening level and boundary conditions. The friction coefficient is updated according to the latest water-bearing state zoning and filling level and the representative value of the unfavorable side is taken. If the behavior consistency value range meets the preset threshold, the final regional geological evaluation model is output. The model includes the final fracture grouping results, fault strong correlation zone, corrected stress distribution and deformation modulus, coordination parameter fusion weight, deformation trend calibration results and spatial location of key control sections.

[0071] like Figure 2 As shown, a subsurface structure intelligent analysis system based on multi-source geological data fusion is also provided to implement the steps of the aforementioned subsurface structure intelligent analysis method based on multi-source geological data fusion. The system includes: The data acquisition and grouping module collects raw geological data and applies a hierarchical clustering algorithm to group scale fractures and interlayer faults to obtain a preliminary stiffness distribution map of scale rock mass units. The deformation mode analysis module uses the finite element analysis method to simulate the stress transfer path between elements based on the preliminary stiffness distribution map of the scale rock mass unit, and determines the deformation mode and failure path of the scale geological block. The iterative optimization module, if the deformation pattern of the scale geological block deviates from the preset threshold by more than a specified range, will reapply the hierarchical clustering algorithm by iteratively adjusting the grouping parameters to obtain an optimized scale rock mass unit stiffness distribution map. The local mechanical response analysis module obtains the optimized scale rock mass element stiffness distribution map, and then uses the finite element analysis method to calculate the local stress concentration and uneven settlement distribution for key parts, including the weak interlayer of the dam foundation, and to determine the location of potential sliding surfaces. The overall mechanical response modeling module, when determining the location of potential sliding surfaces, integrates the deformation patterns of scale geological blocks with local stress concentration data, and applies a multi-scale nested simulation algorithm to construct an overall mechanical response model, thereby obtaining simplified prediction results of the macroscopic behavior of underground structures. The consistency verification module extracts stiffness gradient and slip compatibility from the simplified macroscopic behavior prediction results of geological structures, and uses finite element analysis to verify the matching degree with the actual monitoring data to determine the consistency of the model. The feedback evaluation module, for the determined model coordination consistency, if the matching degree is lower than the preset threshold, will regroup the scale fractures and interlayer faults through feedback loop to obtain the final regional geological evaluation model.

[0072] In the above implementation, the data acquisition and grouping processing module receives raw geological data from the field exploration and monitoring system, completes coordinate unification, anomaly removal, and missing data completion, and establishes a multi-feature input set according to the geometric distribution of fractures and interlayer fault characteristics. Subsequently, hierarchical clustering is performed to obtain fracture groups and fault influence zones. The fracture grouping results and fault influence areas are further integrated to form scale rock mass unit boundaries, and material stiffness and related mechanical parameters are extracted within these boundaries to form a preliminary stiffness distribution map. The output of this module serves as the basis for subsequent material field and structural zoning calculations of the entire system.

[0073] The deformation mode analysis module takes a preliminary stiffness distribution map and element boundaries as input, constructs a finite element mesh, generates inter-element connection descriptions, and solves for stress transfer paths and displacement responses after loading geological load data. Based on the main stress transfer channels, interfaces with significant differential displacements, and strain exceedance sequences, it outputs the deformation modes and failure paths of the geological blocks at the scale. The output of this module includes not only the mode type but also spatially connected links and key control segments, providing locatable feedback signals for iterative optimization and local analysis.

[0074] The iterative optimization module receives the output from the deformation mode analysis module and compares it with a preset threshold system to generate deviation detection results. These deviation detection results drive iterative adjustments to the grouping parameters. During the iteration process, feature weights, truncation levels, and fault boundary constraint strengths are directionally corrected. After each iteration, hierarchical clustering and rapid mechanical verification are restarted until the deviation converges or the stability criterion is met. This module outputs an optimized scale rock mass element stiffness distribution map, enabling subsequent calculations to be based on a material field that more closely reflects the actual response.

[0075] The local mechanical response analysis module takes the optimized stiffness distribution map as input to identify the weak interlayers in the dam foundation and their neighborhood. It extracts material and interface property parameters to construct a local geological model, generates a local finite element mesh, and applies boundary conditions and load distributions consistent with the global model. It introduces hydraulic seepage parameters to calculate the stress field changes after pore pressure involvement, extracts local stress concentration values, and couples them with uneven settlement distribution for discrimination. Parameter correction is triggered by a settlement deviation threshold, and the potential sliding surface location is determined by combining failure path simulation under settlement equilibrium conditions. The output of this module provides key constraints for the overall modeling stage.

[0076] The overall mechanical response modeling module receives deformation mode and local stress concentration data, constructs a multi-scale nested overall mechanical response model, and achieves closed-loop coupling through boundary transfer and equivalent backpropagation between the global model and local sub-models. It outputs simplified macroscopic behavior prediction results for underground structures, and incorporates groundwater level influence parameters into the macroscopic results to form the preliminary location and correction results of the sliding surface. This module's output provides a unified macroscopic prediction caliber for consistency verification.

[0077] The consistency verification module extracts the stiffness gradient distribution and slip compatibility relationship from the macroscopic behavior prediction results, constructs a description of the mechanical response law, and forms a finite element comparison input dataset under the constraints of geological environmental variables. After solving for the preliminary comparison results, it calculates the deviation value by aligning with the actual monitoring data. By adjusting the slip compatibility parameters and the stiffness gradient weight ratio, it obtains the distribution range of behavior consistency values ​​and determines whether the model has achieved compatibility consistency based on a preset threshold. This module outputs the compatibility consistency conclusion and prediction adjustment scheme items, providing trigger conditions and correction directions for the feedback evaluation module.

[0078] The feedback evaluation module receives the consistency conclusions. When the matching degree falls below a preset threshold, a feedback loop is initiated. This loop, centered on fracture grouping iteration and fault reconstruction correlation, updates stress field data and interface friction parameters in tandem. The updated fracture groups, fault strong correlation zones, corrected stress distribution, and deformation modulus are then fed back to the preceding modules for end-to-end recalculation. This process continues until the behavioral consistency value reaches the threshold requirement, at which point the final regional geological evaluation model is output. This module enables the system to adaptively correct itself, ensuring that the final output has a verifiable consistency basis.

[0079] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for intelligent analysis of underground structures based on multi-source geological data fusion, characterized in that, include: S1. By collecting raw geological data and applying a hierarchical clustering algorithm to group scale fractures and interlayer faults, a preliminary stiffness distribution map of scale rock mass units is obtained. S2. Based on the preliminary stiffness distribution map of the scale rock mass unit, the finite element analysis method is used to simulate the stress transfer path between units and determine the deformation mode and failure path of the scale geological block. S3. If the deformation pattern of the scale geological block deviates from the preset threshold by more than a specified range, the hierarchical clustering algorithm is reapplied by iteratively adjusting the grouping parameters to obtain an optimized scale rock mass unit stiffness distribution map. S4. After obtaining the optimized scale rock mass unit stiffness distribution map, for key parts including the weak interlayer of the dam foundation, the finite element analysis method is used to calculate the local stress concentration and uneven settlement distribution, and to determine the location of potential sliding surfaces. S5. When determining the location of potential sliding surfaces, by integrating the deformation patterns of scale geological blocks and local stress concentration data, a multi-scale nested simulation algorithm is applied to construct an overall mechanical response model, and a simplified prediction result of the macroscopic behavior of underground structures is obtained. S6. Extract the stiffness gradient and slip compatibility relationship from the simplified macroscopic behavior prediction results of geological structures, and use the finite element analysis method to verify the matching degree with the actual monitoring data to determine the consistency of the model.

2. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S1 includes: Raw geological data were obtained from field exploration. A hierarchical clustering algorithm was applied to the geological data to process the fractures into groups, and fracture grouping results were obtained. Based on the fracture grouping results, inter-layer fault analysis is performed in conjunction with inter-layer fault data, and the fault-affected area is determined by comparing displacement differences; By integrating the results of fault influence areas and fracture grouping, rock mass units of different scales are divided, and unit boundary descriptions are obtained. Obtain the stiffness parameters within the element boundary description, map the initial stiffness distribution, and generate a draft distribution map; Stability indices were extracted from the draft distribution map, and rock mass units were grouped and optimized to obtain a preliminary stiffness distribution map of the rock mass units at different scales.

3. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S2 includes: Scale rock mass element parameters are obtained from the preliminary stiffness distribution map, and the mesh structure is generated using the finite element method to obtain the inter-element connection description. For the inter-unit connection description, input geological load data to simulate stress transfer and determine path distribution characteristics; The block deformation response is calculated based on path distribution characteristics to determine the mode type; Obtain the strain threshold under the mode type, perform strain interaction analysis on the rock mass element, and obtain the failure path sequence; By integrating block stability indices based on the failure path sequence, the deformation patterns and failure paths of scale geological blocks can be determined.

4. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S3 includes: The threshold deviation detection result is obtained by comparing the deformation mode with the preset threshold. The deviation data is iteratively adjusted by grouping parameters to obtain the adjusted parameter set. For the adjusted parameter set, a hierarchical clustering algorithm is used to restart the grouping process of the rock mass unit data to generate preliminary rock mass unit groups; Based on the preliminary rock mass unit grouping, the strain threshold is obtained from the geological load data to determine the stress balance adjustment scheme between units; By adjusting the stress balance between elements, the load data is processed using simulation input to determine changes in the mesh connectivity description; By obtaining the changes in mesh connectivity description and integrating the failure path sequence and stability index through weighted calculation, an optimized scale rock mass element stiffness distribution map is obtained.

5. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S4 includes: By optimizing the stiffness distribution map of the rock mass unit at the scale, the material property parameters of the weak interlayer area of ​​the dam foundation are obtained, and local geological model data is obtained. For local geological model data, a mesh model is constructed using the finite element analysis method to determine boundary conditions and load distribution; The hydraulic seepage influence parameters are input into the grid model. These parameters are obtained in advance from the geological exploration data of the dam foundation. The stress field changes in the weak interlayer of the dam foundation are calculated. The calculation process uses the stress balance equation to integrate the seepage pressure distribution and obtain the local stress concentration value. Based on the local stress concentration value, integrate the uneven settlement distribution data. If the settlement deviation exceeds the preset threshold, adjust the model parameters and determine the settlement equilibrium state. By obtaining the settlement equilibrium state and combining it with the failure path simulation, the location of the potential sliding surface can be determined.

6. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S5 includes: By using the deformation model of a geological block at scale, local stress concentration data can be obtained, and relevant parameters of the deformation model can be integrated to determine the distribution of the initial response of the underground structure. For the initial response distribution, a multi-scale nested simulation algorithm is adopted. By integrating deformation mode-related parameters and local stress data, an overall mechanical response model is constructed to obtain a simplified macroscopic behavior distribution. Geological block information is obtained from the macroscopic behavior distribution, groundwater level influence parameters are calculated, and the preliminary location of potential sliding surfaces is determined. Obtain the preliminary position, adjust the relevant parameters of the deformation mode based on the local stress concentration data, and determine the correction value of the sliding surface position; Based on the corrected values, relevant data on the distribution of macroscopic behavior are integrated to obtain the prediction results of the macroscopic behavior of underground structures.

7. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that: S6 includes: From the geological structure prediction results, relevant data on stiffness gradient distribution are obtained. Based on the stiffness gradient distribution and the changing characteristics of underground deformation trends, a preliminary description of mechanical response laws is constructed, and the initial distribution characteristics of mechanical response laws are obtained. Based on the initial distribution characteristics of the mechanical response law, the specific value range of the slip compatibility parameter is obtained. By mapping the slip compatibility parameter with the local stress distribution, the boundary conditions of the mechanical response law are adjusted, and the corrected response distribution range is determined. Based on the corrected response distribution range, relevant information on macroscopic behavioral characteristics is integrated, and an input dataset for the finite element analysis method is constructed under the influence of geological environmental variables to obtain preliminary comparative results of the finite element analysis method. By comparing the preliminary results, the deviation value of the monitoring data comparison is obtained. Based on the requirements of data matching accuracy, the weight ratio of the slip coordination parameter and the stiffness gradient distribution is adjusted to determine the distribution range of the behavior consistency value. If the distribution range of behavioral consistency values ​​meets the preset threshold, then by combining the relevant data of local stress distribution and underground deformation trend, the final geological structure prediction and adjustment scheme is generated, and the verification conclusion of macroscopic behavioral characteristics is determined.

8. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 1, characterized in that, It also includes S7, which addresses the consistency of a determined model. If the matching degree is lower than a preset threshold, a feedback loop is used to regroup scale fractures and interlayer faults to obtain the final regional geological evaluation model, specifically including: Specific data on fault reconstruction correlation are obtained from the fracture grouping iteration. The fracture grouping iteration performs multiple cyclic classifications based on the geometric distribution and stress field data of scale fractures. The fault reconstruction correlation is calculated using the inter-layer displacement vector and friction coefficient. Based on the results of the matching degree threshold verification, an initial framework for feedback loop optimization is constructed to obtain the preliminary distribution of model consistency adjustment. For the preliminary distribution, relevant information from the scale feature simulation is integrated. The scale feature simulation uses the finite difference method to process fracture density and permeability data. Through the interlayer interaction correction process, the interaction vector is adjusted according to the displacement gradient and shear modulus. The boundary conditions of the deviation value correction process are adjusted to determine the correction range of the integrated geological evaluation.

9. The intelligent analysis method for underground structures based on multi-source geological data fusion according to claim 8, characterized in that: The S7 also includes: Based on the corrected value range, the input dataset for constructing the regional model is obtained. The input dataset contains the corrected stress distribution and deformation modulus. Under the influence of the coordination parameter fusion, where the coordination parameter fusion combines the stiffness gradient and slip parameter through a weighted average method, the deformation trend calibration operation is performed to obtain the final scheme optimization comparison results. By comparing the results, the range of behavioral consistency values ​​is determined. If it is lower than the preset threshold, the fracture grouping iteration and fault reconstruction association are re-executed. The re-execution includes updating the stress field data and friction coefficient to obtain the final regional geological evaluation model.

10. An intelligent analysis system for underground structures based on multi-source geological data fusion, used to implement the steps of the intelligent analysis method for underground structures based on multi-source geological data fusion as described in any one of claims 1-9, characterized in that, The system includes: The data acquisition and grouping module collects raw geological data and applies a hierarchical clustering algorithm to group scale fractures and interlayer faults to obtain a preliminary stiffness distribution map of scale rock mass units. The deformation mode analysis module uses the finite element analysis method to simulate the stress transfer path between elements based on the preliminary stiffness distribution map of the scale rock mass unit, and determines the deformation mode and failure path of the scale geological block. The iterative optimization module, if the deformation pattern of the scale geological block deviates from the preset threshold by more than a specified range, will reapply the hierarchical clustering algorithm by iteratively adjusting the grouping parameters to obtain an optimized scale rock mass unit stiffness distribution map. The local mechanical response analysis module obtains the optimized scale rock mass element stiffness distribution map, and then uses the finite element analysis method to calculate the local stress concentration and uneven settlement distribution for key parts, including the weak interlayer of the dam foundation, and to determine the location of potential sliding surfaces. The overall mechanical response modeling module, when determining the location of potential sliding surfaces, integrates the deformation patterns of scale geological blocks with local stress concentration data, and applies a multi-scale nested simulation algorithm to construct an overall mechanical response model, thereby obtaining simplified prediction results of the macroscopic behavior of underground structures. The consistency verification module extracts stiffness gradient and slip compatibility from the simplified macroscopic behavior prediction results of geological structures, and uses finite element analysis to verify the matching degree with the actual monitoring data to determine the consistency of the model. The feedback evaluation module, for the determined model coordination consistency, if the matching degree is lower than the preset threshold, will regroup the scale fractures and interlayer faults through feedback loop to obtain the final regional geological evaluation model.

Images