Multi-medium three-dimensional modeling and geotechnical mechanics response coupling method

Abstract

A multi-medium three-dimensional modeling and geotechnical mechanics response coupling method, comprising: establishing a collapse risk assessment model according to the time evolution data of the deformation accumulation, if the local deformation rate exceeds the preset threshold, it is judged as a high-risk area, if the deformation accumulation reaches the critical value, the potential collapse position is determined; the space prediction model is constructed through the high-risk area and the potential collapse position information, the risk assessment algorithm is used to quantitatively calculate the collapse probability of different areas, and the spatial distribution prediction result of karst collapse is obtained; the time and space information output by the early warning mechanism is used to update the numerical model parameters, the prediction accuracy is optimized and adjusted through the feedback correction mechanism, and the karst collapse prediction technology system suitable for complex geological conditions is obtained. The present application realizes the accurate prediction and early warning of karst collapse through the mode of data driving and model coupling, provides technical support for the prevention and control of karst collapse under complex geological conditions, and has important practical value.

Description

Technical Field

[0001] This invention relates to the field of information technology, and in particular to a multi-media three-dimensional modeling and geotechnical mechanical response coupling method. Background Technology

[0002] Karst collapse, as one of the most destructive geological hazards, is characterized by its suddenness and high degree of harm, making accurate prediction and effective prevention and control extremely urgent.

[0003] Current research on karst collapse often relies on single-discipline theories for analysis, lacking comprehensive research methods that integrate multiple disciplines. Existing prediction methods mainly depend on empirical judgment and qualitative analysis, making it difficult to achieve accurate quantitative assessments, and their applicability and reliability under complex geological conditions are significantly insufficient.

[0004] The complexity of karst collapse mechanisms stems from the strong coupling effect between groundwater dynamics and the mechanical response of soil and rock. Fluctuations in groundwater level not only alter the stress state of the soil and rock mass but also trigger a redistribution of pore water pressure and effective stress. This hydraulic-mechanical coupling process directly controls the deformation characteristics of the soil. The positive pressure exerted by rising water levels and the vacuum erosion effect created by falling water levels generate complex stress concentration phenomena around karst caves and conduits near the soil-rock interface. This stress concentration further exacerbates uneven soil deformation, particularly at the top of karst caves and on both sides of karst conduits, where the accumulation of localized deformation may ultimately trigger a collapse event.

[0005] How to establish a numerical model that can accurately describe the coupling mechanism between groundwater dynamics and soil mechanics response, realize the quantitative simulation of the mechanical effects of soil-rock interface under different groundwater level changes, and reveal the intrinsic mechanism of karst collapse through accurate calculation of effective stress and soil deformation, thereby providing reliable technical support for the spatial prediction of karst collapse, has become a key issue that urgently needs to be solved in the field of karst geological disaster prevention and control. Summary of the Invention

[0006] This invention provides a multi-media three-dimensional modeling and geotechnical mechanical response coupling method, mainly including:

[0007] Data on groundwater level monitoring and geological exploration in karst areas are acquired. The monitoring data is then processed by a data preprocessing module to perform time-series standardization and spatial interpolation, resulting in a continuous spatiotemporal distribution of water level change sequences and geometric parameter information of the soil-rock interface.

[0008] A numerical model of groundwater seepage field was established based on the water level change sequence. The fluid-structure interaction algorithm was used to calculate the pore water pressure distribution and seepage force field under different water level conditions, and the hydraulic gradient and seepage velocity distribution characteristics near the soil-rock interface were obtained.

[0009] A stress field analysis model for soil and rock mass is constructed using the pore water pressure distribution data. If the water level rises, the top support stress effect is calculated; if the water level falls, the vacuum absorption stress effect is calculated, thus determining the stress concentration areas around the karst cave and on both sides of the karst conduit.

[0010] The stress distribution data of the stress concentration area is used to calculate the soil deformation. The displacement field and strain field of the soil-rock interface are solved by the finite element method to obtain the spatial distribution and temporal evolution law of local deformation accumulation.

[0011] A collapse risk assessment model is established based on the time-series evolution data of the accumulated deformation. If the local deformation rate exceeds a preset threshold, it is judged as a high-risk area. If the accumulated deformation reaches a critical value, the potential collapse location is determined.

[0012] A spatial prediction model is constructed using the high-risk areas and potential collapse location information. A risk assessment algorithm is used to quantitatively calculate the collapse probability of different areas, and the spatial distribution prediction results of karst collapse are obtained.

[0013] A dynamic early warning system is established based on the spatial distribution prediction results and real-time monitoring data. If the trend of the monitoring parameters changes is consistent with the output of the prediction model, the early warning mechanism is triggered to determine the time window and scope of impact of the collapse.

[0014] The spatiotemporal information output by the aforementioned early warning mechanism is used to update the parameters of the numerical model, and the prediction accuracy is optimized and adjusted through a feedback correction mechanism to obtain a karst collapse prediction technology system adapted to complex geological conditions.

[0015] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects:

[0016] This invention acquires groundwater level monitoring data and geological exploration data to establish a numerical model of the groundwater seepage field and a stress field analysis model of the soil and rock mass. It calculates the pore water pressure distribution and seepage force field under different water level conditions, identifies stress concentration areas, and calculates soil deformation to obtain the spatial distribution and temporal evolution of local deformation accumulation. Based on this, a collapse risk assessment model is established to spatially predict and probabilistically calculate high-risk areas and potential collapse locations, constructing a dynamic early warning system. This invention, through data-driven and model-coupled approaches, achieves accurate prediction and early warning of karst collapse, providing technical support for the prevention and control of karst collapse under complex geological conditions, and has significant practical value. Attached Figure Description

[0017] Figure 1 This is a flowchart of a multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to the present invention.

[0018] Figure 2This is a schematic diagram of a multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to the present invention. Detailed Implementation

[0019] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this specification, and not all embodiments. Based on the embodiments in this specification, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this specification.

[0020] like Figure 1-2 This embodiment of a multi-media three-dimensional modeling and geotechnical mechanical response coupling method specifically includes:

[0021] Step S101: Obtain groundwater level monitoring data and geological exploration data in karst areas. Perform time-series standardization and spatial interpolation on the monitoring data through the data preprocessing module to obtain a continuous spatiotemporal distribution of water level change sequence and geometric parameter information of soil-rock interface.

[0022] Groundwater level monitoring data and geological exploration data were acquired from karst areas using a data acquisition system to obtain a raw dataset. A data preprocessing module was used to perform time-series standardization on the groundwater level monitoring data in the raw dataset, generating time-series consistent water level change data. For the water level change data, spatial interpolation methods were used to fill in missing areas, resulting in a continuous spatiotemporally distributed water level change sequence. Based on the geological exploration data, the distribution characteristics of the soil-rock interface were extracted, and the geometric parameters of the soil-rock interface were determined. If the water level change sequence showed abnormal fluctuations, the geometric parameters of the soil-rock interface were used to determine whether there were areas at risk of karst cave collapse. After obtaining the judgment results, the spatial distribution pattern of stratigraphic unconformities was analyzed based on the distribution of carbonate rocks and fracture development characteristics within the risk area to determine the potential collapse impact range. By analyzing the water level change trend and geological characteristics within the potential collapse impact range, a distribution layer of karst cave collapse risk was generated.

[0023] Specifically, in the acquisition and processing of groundwater level monitoring data and geological exploration data in karst areas, groundwater level data was first obtained from 50 monitoring wells distributed throughout the study area using an automated data acquisition system. The data spanned from January to December 2022, with water level values ​​recorded daily. For example, a monitoring well recorded a water level of 125.3 meters on June 1, 2022. Next, borehole data at the soil-rock interface was extracted using a geological exploration database, collecting information from 100 borehole locations. Each location included depth and lithological descriptions; for example, a location with a depth of 30.5 meters and a lithology of limestone was identified. Following this, the data preprocessing stage involved time-series standardization of the water level monitoring data. The Z-score standardization method was used to convert the water level data from each monitoring well into a standardized sequence with a mean of 0 and a standard deviation of 1. The calculation formula is as follows: Where X is the original water level value, The mean of the time series. The standard deviation is calculated as follows: for example, the water level of a well on June 1st was 125.3 meters, with a mean of 120.0 meters and a standard deviation of 5.0 meters, resulting in a standardized value of 1.06. Simultaneously, spatial interpolation was performed on the spatially unevenly distributed water level data. The Kriging interpolation algorithm was used to generate a 10m × 10m resolution water level grid within the study area using ArcGIS software. The interpolation results showed that the water level at the center point of a certain area was 123.8 meters, with an error range controlled within ±0.5 meters, forming a continuous spatiotemporal sequence of water level changes. For the extraction of geometric parameters at the soil-rock interface, geological exploration borehole data was used, combined with triangular mesh interpolation (TIN) to construct a three-dimensional interface model. The calculated average interface depth was 28.7 meters, and the maximum slope was 15.2 degrees. The model accuracy was evaluated through cross-validation, with the mean error controlled within 0.3 meters. These processing results provide a reliable basis for subsequent hydrogeological analysis. At the same time, correlation analysis with regional rainfall data (e.g., an average annual rainfall of 800 mm) revealed a positive correlation between water level changes and rainfall, with a correlation coefficient of 0.75. This further validates the rationality of the data processing and forms a complete logical chain from data acquisition to analysis.

[0024] Step S102: Establish a numerical model of groundwater seepage field based on the water level change sequence, and use fluid-structure interaction algorithm to calculate the pore water pressure distribution and seepage force field under different water level conditions to obtain the hydraulic gradient and seepage velocity distribution characteristics near the soil-rock interface.

[0025] Based on water level change sequences and geological exploration data, a three-dimensional numerical model of the groundwater seepage field was constructed. The finite element method was used to calculate the pore water pressure distribution, obtaining the seepage force field characteristics. Using these seepage force field characteristics, combined with the geometric parameters of the soil-rock interface, the hydraulic gradient and seepage velocity distribution were calculated to determine the seepage path within the overlying carbonate rock area. If the seepage velocity exceeds a preset threshold, the unconformity surface and fracture development characteristics were extracted to generate distribution data for high-risk seepage areas. Based on the distribution data of high-risk seepage areas, the spatial correlation between carbonate rock distribution and fracture networks was analyzed to obtain the boundary information of potential impact areas for cavern collapse. Using the boundary information of potential impact areas, combined with the hydraulic gradient distribution, a spatial distribution layer of collapse risk was generated using Kriging interpolation. Based on the spatial distribution layer of collapse risk, the superimposed areas of fracture development characteristics and unconformity surfaces were extracted to determine priority monitoring areas. Using monitoring data from these priority areas, the seepage field model parameters were updated to generate dynamically adjusted collapse risk distribution data.

[0026] Specifically, based on the groundwater level change sequence in karst areas, a groundwater seepage field model was established using the numerical simulation software MODFLOW. The study area was set as a 5 km × 5 km rectangular grid with a resolution of 50 m × 50 m and a total of 10,000 grid cells. The model input was groundwater level data from January to December 2023, with a time step of 1 day. For example, the groundwater level at a monitoring point was 128.2 meters on March 15, 2023. A fluid-structure interaction algorithm was employed, using the finite difference method to solve the groundwater flow equation, coupled with soil deformation effects, to calculate the pore water pressure distribution. In the algorithm, the pore water pressure P was calculated using Darcy's law, with the formula: ,in With a water density of 1000 kg / m³, Given the gravitational acceleration of 9.8 m / s², h as the hydraulic head, and z as the elevation, for example, at a point with a hydraulic head of 130.5 m and an elevation of 20 m, the calculated value is P = 1087800 Pascals. The seepage force field is generated through the hydraulic gradient. Calculations were performed, and the seepage velocity was obtained by combining the permeability coefficient k (0.001 m / s for limestone). .

[0027] For example, at a certain grid point, the hydraulic gradient is 0.02, and the seepage velocity is 0.00002 m / s. For the soil-rock interface, based on an existing 3D geological model, the hydraulic gradient and seepage velocity of the grid near the interface are extracted. Analysis shows that the average hydraulic gradient at the interface is 0.015, and the maximum seepage velocity is 0.00003 m / s. The model boundary conditions are set as fixed head boundaries, with a head of 135 meters at the northern boundary and 120 meters at the southern boundary. The calculation results are post-processed using a Python script to generate a seepage velocity distribution map, showing a higher velocity in the central region, consistent with the distribution of geological structural fractures. By comparing with measured water levels, the model error is controlled within ±0.4 meters, verifying the model's reliability. Further correlation with the regional permeability coefficient distribution reveals that the high-permeability zone (k=0.002 m / s) overlaps with the high-seepage-velocity zone, with a correlation coefficient of 0.82, revealing the coupling relationship between the seepage field and geological conditions, providing data support for subsequent water resource management.

[0028] Step S103: Construct a stress field analysis model for soil and rock mass using the pore water pressure distribution data. If the water level rises, calculate the top support stress effect; if the water level falls, calculate the vacuum absorption stress effect to determine the stress concentration areas around the karst cave and on both sides of the karst conduit.

[0029] Based on pore water pressure distribution data, a three-dimensional stress field model of the rock and soil mass is constructed using the finite element method. The stress distribution characteristics under different water level changes are calculated, yielding spatial distribution data of stress concentration areas. Using this spatial distribution data, combined with the distribution characteristics of stratigraphic unconformities and fracture surfaces, spatial overlay analysis is employed to determine the overlapping areas between stress concentration areas and overlying carbonate rock distributions, thus obtaining boundary information for high-risk areas. If the boundary information of high-risk areas overlaps with karst conduit distributions, geostatistical analysis is used to calculate the stress gradient distribution within the stress concentration areas, identifying areas of stress gradient anomalies. Based on these anomaly areas and the geological features surrounding the karst caves, Kriging interpolation is used to generate a spatial distribution layer of the stress gradient anomaly areas, identifying potential collapse risk areas. Through these potential collapse risk areas, the overlay characteristics of fracture surface distribution and stratigraphic unconformities are extracted. Spatial clustering analysis is used to obtain a priority ranking of high-risk areas, determining priority monitoring areas. Based on the stress distribution characteristics of the priority monitoring areas, time series analysis is used to update the stress field model parameters, generating dynamically adjusted stress distribution data. By using dynamically adjusted stress distribution data, combined with the characteristics of top support stress and vacuum erosion stress, the finite difference method is used to calculate the stress change trend in high-risk areas, thus obtaining the dynamic distribution of potential collapse areas.

[0030] Specifically, based on pore water pressure distribution data in karst areas from January to December 2024, a stress field model of the soil and rock mass was constructed using the finite element analysis software ABAQUS. The study area was set as a 4 km × 4 km rectangular grid with a resolution of 40 m × 40 m, and a total of 10,000 grid cells. The input data was the daily pore water pressure; for example, the pore water pressure at a monitoring point on May 10, 2024, was 1,050,000 Pascals. First, the stress field was calculated using the pore water pressure data, employing the Mohr-Coulomb constitutive model, with the following formula: ,in For effective stress, Let the initial stress be 500,000 Pascals. is the pore pressure coefficient (taken as 0.8), and P is the pore water pressure.

[0031] For example, at a point P = 1050000 Pascals, the calculation is... = 1,340,000 Pascals. When the water level rises by 5 meters (the water head increases to 132 meters), the stress is calculated using the formula for top support stress. ( =1000 kg / m³ =9.8 m / s², =5 meters), calculate the stress increment of the top support. =49000 Pascals, superimposed on the stress field. When the water level drops by 3 meters (water head drops to 124 meters), the vacuum absorption stress formula is used. (k is the absorption coefficient, taken as 0.5), calculate the absorption stress. =-14700 Pascal, stress distribution updated. Mesh stress data was extracted around the cave (radius 30 m) and on both sides of the karst conduit (width 20 m). Analysis revealed stress concentration at the top of the cave, with a maximum stress of 1450000 Pascals, and a stress gradient of 5000 Pascals / m on both sides of the conduit. The stress field data was processed using a Python script to generate a stress distribution cloud map, showing a correlation between the stress concentration area around the cave and the location of geological fissures, with a correlation coefficient of 0.78. Boundary conditions were set to a fixed stress of 1000000 Pascals in the north and 800000 Pascals in the south to ensure model convergence. The analysis results indicate a high correlation between stress concentration areas and the degree of karst development, providing a basis for karst collapse risk assessment.

[0032] Step S104: The stress distribution data of the stress concentration area is used to calculate the soil deformation. The displacement field and strain field of the soil-rock interface are solved by the finite element method to obtain the spatial distribution and temporal evolution law of local deformation accumulation.

[0033] Using the finite element method (FEM) based on stress distribution data, the displacement and strain fields at the soil-rock interface are calculated to obtain the spatial distribution characteristics and temporal evolution of local deformation accumulation. Based on the spatial distribution characteristics of local deformation accumulation, and combined with the planar distribution of karst cave collapse areas within the overburdened carbonate rock distribution area, spatial analysis methods are used to extract the superimposed characteristics of stratigraphic unconformities, fracture surfaces, and fracture zone development to identify potential high-risk deformation areas. Using these high-risk areas, and combining infill material distribution data, Kriging interpolation is used to generate a spatial distribution layer of deformation gradients within the high-risk areas, identifying areas with abnormal deformation gradients. If an abnormal deformation gradient area overlaps with a karst cave collapse area, geostatistical analysis is used to calculate the deformation rate distribution within the abnormal area, obtaining an abnormal deformation rate area. Based on the abnormal deformation rate areas, and combining the distribution characteristics of stratigraphic unconformities and fracture surfaces, spatial clustering analysis is used to determine the priority ranking of high-risk areas, identifying priority monitoring areas. Using the deformation rate distribution of the priority monitoring areas, time series analysis is used to update the displacement and strain field model parameters, generating dynamically adjusted deformation distribution data. Based on the dynamically adjusted deformation distribution data, combined with the distribution of filling material and the development characteristics of the fracture zone, the finite difference method is used to calculate the deformation change trend of high-risk areas and determine the dynamic distribution of potential collapse areas.

[0034] Specifically, based on stress distribution data of karst areas from June to December 2024, soil deformation was calculated using the finite element analysis software COMSOL Multiphysics. Displacement and strain field models of the soil-rock interface were constructed. The study area was set as a 5 km × 5 km square grid with a resolution of 50 m × 50 m, and a total of 10,000 grid cells. The input data was the daily stress distribution; for example, the stress value at a monitoring point on August 15, 2024, was 1,200,000 Pascals. A linear elastic constitutive model was adopted, with the following formula: ,in For strain, E is the elastic modulus of the soil (taken as 2.5 × 10^9 Pascals). This is the input stress.

[0035] For example, a certain point =1200000 Pascals, calculate the strain as follows: =0.00048. The displacement field was solved using the finite element method, employing a quadrilateral mesh with mesh elements satisfying a minimum angle condition (greater than 30 degrees). Boundary conditions were set as follows: a fixed displacement of 0 meters in the west and a uniform load of 500,000 Pascals applied in the east. The maximum displacement at a certain point was calculated to be 0.015 meters, located near the soil-rock interface. By iteratively solving the displacement field and combining it with a time step (updated daily), the temporal evolution was analyzed, revealing that the cumulative displacement increased by 0.008 meters from September to November, correlated with rainfall data (average monthly rainfall of 200 mm), with a correlation coefficient of 0.82. The strain field analysis used strain gradient calculation, with the formula: Using Python scripts to process grid data, strain concentration areas at the soil-rock interface (10 meters thick) were extracted. The maximum strain of 0.0006 was found to occur at the interface center, with an elliptical spatial distribution consistent with the geological fault strike (angle deviation less than 5 degrees). Displacement and strain distribution contour maps generated using MATLAB showed that local deformation was concentrated in the area 20 meters below the soil-rock interface, consistent with the distribution of karst fissures (consistency 0.85). To ensure model convergence, the Newton-Raphson iteration method was used, with residuals controlled within 10^-6. The analysis results provide data support for predicting surface subsidence in karst areas.

[0036] Step S105: Establish a collapse risk assessment model based on the time-series evolution data of the accumulated deformation. If the local deformation rate exceeds a preset threshold, it is judged as a high-risk area. If the accumulated deformation reaches a critical value, the potential collapse location is determined.

[0037] By using time-series evolution data and employing time-series decomposition methods, the periodic and trend components of local deformation rates are extracted to obtain the local deformation rate distribution. If the trend component of the local deformation rate exceeds a preset threshold, a spatial distribution layer of high-risk areas is generated using spatial interpolation methods to determine the scope of high-risk areas. For the high-risk area scope, combined with deformation accumulation data, a threshold comparison method is used. If the accumulation exceeds a critical value, a spatial positioning algorithm is used to determine the set of potential collapse locations. Based on the set of potential collapse locations and combined with carbonate rock distribution area data, spatial overlay analysis is used to extract the overlapping areas of stratigraphic unconformities and fracture surfaces, obtaining geologically vulnerable areas. Through the geologically vulnerable areas, combined with fracture zone characteristic data, density analysis methods are used to calculate the fracture zone distribution density within the vulnerable areas, identifying high-density fracture zone areas. For high-density fracture zone areas, combined with infill material type data, a classification algorithm is used to extract the distribution ratios of sand, clay, and gravel, obtaining the infill material distribution characteristics. Based on the distribution characteristics of the filling material, a weighted overlay analysis was adopted to integrate data from areas with fragile geological structures and high-density fracture zones, generating a priority distribution of collapse risk and determining priority monitoring areas.

[0038] Specifically, based on surface deformation monitoring data of karst areas from June to December 2024, a collapse risk assessment model was constructed using the Python programming language. The analysis area was set as a 4 km × 4 km square region with a resolution of 40 m × 40 m and a total grid size of 10,000. Daily deformation rate data were acquired using remote sensing InSAR technology; for example, the deformation rate at a certain monitoring point on July 20, 2024, was 0.002 m / day. First, a time series analysis algorithm (ARIMA model) was used to process the deformation rate data, with model parameters set to (p=2, d=1, q=1). The deformation rate trend for the next 7 days was predicted using least squares fitting; for example, the deformation rate at a certain point on July 27, 2024, was predicted to be 0.0023 m / day. Next, a deformation rate threshold of 0.0025 m / day was set. If a grid point exceeded this threshold for three consecutive days, it was marked as a high-risk area. The calculation results showed that approximately 1200 grid points were marked, accounting for 12% of the total number of grid points. Subsequently, the cumulative deformation was calculated using the following formula: Where D is the cumulative deformation, This represents the daily deformation rate. The time step is 1 day; for example, the cumulative deformation at a point from June 1st to December 31st is 0.42 meters. A critical deformation value of 0.45 meters is set. If the cumulative deformation at a point exceeds this value, it is identified as a potential collapse location. Analysis revealed that approximately 300 grid points meet the criteria, mainly concentrated in the southeastern part of the region. A spatial interpolation algorithm (Kriging method) is used to generate a spatial distribution map of cumulative deformation. The interpolation radius is set to 100 meters, showing that high-risk areas are strip-shaped and correlated with the distribution of underground karst channels, with a correlation coefficient of 0.78. To improve model accuracy, a Support Vector Machine (SVM) algorithm is used to classify high-risk areas. The feature vectors include deformation rate, cumulative deformation, and terrain slope (average slope 5 degrees), achieving a classification accuracy of 0.90. Finally, the model outputs risk distribution data, stored in GeoTIFF format for easy loading and visualization in GIS software, providing support for collapse early warning in karst areas.

[0039] Step S106: Construct a spatial prediction model using the high-risk area and potential collapse location information, and use a risk assessment algorithm to quantitatively calculate the collapse probability of different areas to obtain the spatial distribution prediction results of karst collapse.

[0040] By using data on high-risk areas and potential collapse locations, a kernel density estimation algorithm is employed to calculate the spatial probability distribution of karst collapse, resulting in a collapse risk prediction layer. Based on this layer, the spatial extent of the overburden carbonate rock distribution area is extracted. Combined with the distribution data of stratigraphic unconformities and fracture surfaces, spatial overlay analysis is used to determine the set of geologically vulnerable areas. If the set of geologically vulnerable areas overlaps with the distribution of fracture zones, spatial correlation analysis is used to calculate the fracture zone density in the overlapping area, identifying high-density fracture zone areas. For these high-density fracture zone areas, infill material distribution data is obtained, and a classification algorithm is used to extract the distribution ratios of sand, clay, and gravel, yielding the infill material characteristic distribution. Based on this infill material characteristic distribution, combined with the set of geologically vulnerable areas, a weighted overlay analysis is used to generate a collapse sensitivity distribution layer. Using this layer, a spatial clustering algorithm is employed to determine the set of monitoring areas with the highest priority. For this set of highest-priority monitoring areas, combined with historical collapse event data, time series analysis is used to predict the probability distribution of future collapses, resulting in a dynamic risk distribution.

[0041] Specifically, based on surface deformation monitoring data of karst areas from January to May 2024, a spatial prediction model was constructed using the Python programming language. The analysis area was set as a 5 km × 5 km square region with a resolution of 50 m × 50 m and a total grid size of 10,000. Daily deformation rate data were obtained using InSAR technology; for example, the deformation rate of a monitoring point on March 15, 2024, was 0.0018 m / day. A random forest algorithm was used to process the deformation rate data, with 100 trees and a maximum depth of 10. The model was optimized through cross-validation to predict the deformation rate trend for the next 10 days; for example, the predicted deformation rate of a point on March 25, 2024, was 0.0021 m / day. A risk probability threshold of 0.75 was set. If the predicted probability of a grid point exceeded this value, it was marked as a high-risk area. The calculation results showed that approximately 1,500 grid points were marked, accounting for 15% of the total grid size. Subsequently, the cumulative deformation probability was calculated using the formula... Where P is the cumulative risk probability. For daily predicted probabilities, The weights are set as follows (based on topographic factors, with an average weight of 0.8). For example, the cumulative probability of a certain point up to May 31st is 0.82. A critical probability value of 0.85 is set; if this value is exceeded, it is identified as a potential collapse location. Analysis revealed that approximately 400 grid points meet the criteria, mainly concentrated in the northwest of the region. A spatial distribution map of risk probability is generated using the inverse distance weighted (IDW) interpolation algorithm, with an interpolation radius of 120 meters. The map shows that high-risk areas are patchy and correlated with groundwater level changes, with a correlation coefficient of 0.82. The gradient boosting tree (GBDT) algorithm is used to classify high-risk areas. The feature vector includes deformation rate, geological structure factor (fault density 0.02 faults / km), and soil moisture (average moisture 30%), with a classification accuracy of 0.92. Finally, the model outputs risk probability distribution data, stored in NetCDF format, for easy subsequent analysis and visualization, providing support for collapse prediction in karst areas.

[0042] Step S107: Establish a dynamic early warning system based on the spatial distribution prediction results and real-time monitoring data. If the trend of monitoring parameter changes is consistent with the output of the prediction model, the early warning mechanism is triggered to determine the time window and scope of impact of the collapse.

[0043] By combining spatial distribution data and prediction results, stratigraphic distribution information within the overburdened carbonate rock area is obtained. Combined with real-time monitoring data, abnormal areas of parameter changes are identified, resulting in a preliminary risk distribution range. Based on this preliminary risk distribution range, relevant data on fracture fracturing are acquired for these abnormal areas. Spatial overlay analysis is used to determine the degree of overlap between fracture fracturing and stratigraphic distribution, identifying a set of high-risk areas. If parameter changes within this high-risk area set exceed a preset threshold, the trend within a time window is extracted using real-time monitoring data to obtain the potential collapse period distribution. Based on this potential collapse period distribution and impact range data, the distribution characteristics of the infill material are obtained. A classification algorithm is used to determine the weight of the infill material's influence on collapse, resulting in a risk-weighted distribution. Using this risk-weighted distribution, updated real-time monitoring data is acquired to meet the needs of dynamic early warning. If the updated data matches the trend of the risk-weighted distribution, an early warning mechanism is triggered, identifying priority warning areas. Based on these priority warning areas, combined with stratigraphic distribution and fracture fracturing data, a spatial clustering algorithm is used to extract key monitoring points within the impact range, resulting in a set of key monitoring targets. By focusing on a set of key monitoring targets, a continuous real-time data stream is obtained. If the continuous data stream shows that parameter changes are consistently abnormal, the predicted value of the time window is updated to determine the final warning time range.

[0044] Specifically, based on the need for subsidence prediction in karst areas, when constructing a dynamic early warning system, the first step is to integrate spatial distribution prediction results with real-time monitoring data. A time series analysis algorithm is then used to extract trends from the monitoring parameters. Assuming the monitoring area covers a 10 km × 10 km range with a resolution of 100 m × 100 m and a total of 10,000 grid points, the real-time monitoring data includes surface subsidence rate and groundwater level fluctuations. For example, on June 1, 2024, the subsidence rate of a certain grid point was 0.0025 m / day, and the groundwater level drop rate was 0.1 m / day. Next, a Long Short-Term Memory (LSTM) network algorithm is used to train the system on the monitoring data from the past 30 days, with a time step of 5 days and 64 hidden layer units, to predict the parameter change trends for the next 7 days. For example, the predicted subsidence rate of a certain point on June 7, 2024, is 0.0027 m / day. Then, the predicted trend is compared with the output of the spatial distribution model, and the deviation value between the two is calculated. A deviation threshold of 0.0003 meters / day is set; if this value is exceeded, the trend is considered consistent, triggering an early warning mechanism. Analysis shows that approximately 800 grid points meet the criteria. Further, using a time window analysis algorithm combined with historical collapse event data, a time window of 5 to 10 days is set to calculate the possible time range of collapses. For example, the predicted window for a high-risk area is June 8th to June 12th, 2024. Simultaneously, using a spatial buffer analysis algorithm, with an impact radius of 200 meters, the number of affected grid points is calculated, resulting in approximately 1200 grid points potentially affected. Finally, the early warning information and impact range data are integrated to automatically generate an early warning report, stored in GeoJSON format, and pushed to the relevant management platform through a system interface, ensuring timely information transmission and forming a complete logical chain from data analysis to early warning triggering.

[0045] Step S108: Update the numerical model parameters using the spatiotemporal information output by the early warning mechanism, and optimize and adjust the prediction accuracy through a feedback correction mechanism to obtain a karst collapse prediction technology system adapted to complex geological conditions.

[0046] The spatiotemporal information obtained through the early warning mechanism is used to dynamically adjust the parameters of the numerical model. Information processing is performed on the data characteristics under complex geological environments to obtain an updated set of model parameters. Based on the updated set of model parameters, a feedback correction mechanism is used to calibrate the prediction accuracy. Combining planar distribution data and carbonate rock regional characteristics, the range of calibrated prediction deviations is determined. If the calibrated prediction deviation range exceeds a preset threshold, the distribution characteristics of potential karst collapse areas are extracted through the fusion processing of spatiotemporal information and geological environment data, resulting in the high-risk area delineation. For the high-risk area delineation, data on fracture distribution and fractured zone areas are obtained. Spatial analysis tools are used to compare the overlap between fracture distribution and planar distribution to determine the set of potential collapse areas. If the fracture distribution characteristics within the potential collapse area set are highly correlated with the fractured zone data, auxiliary analysis is performed using infill material distribution information, combined with geological environment characteristics, to determine the scope of key areas of concern. Based on the scope of key areas of concern, real-time data streams are used to monitor the dynamic changes of karst collapses. Combined with the output of the feedback correction mechanism, an updated risk distribution layer is obtained. By using the updated risk distribution layer and combining it with the spatiotemporal information output of the early warning mechanism, a logistic regression algorithm is used to classify the risk distribution layer and determine the final set of high-risk warning areas.

[0047] Specifically, based on the need for subsidence prediction in karst areas, a prediction technology system adapted to complex geological conditions is constructed. This system optimizes numerical model parameters and improves prediction accuracy by integrating spatiotemporal information output from the early warning mechanism. First, spatiotemporal data generated by the early warning system is collected. For example, the predicted surface subsidence rate for a certain area on June 10, 2024, is 0.0030 m / day, and the groundwater level change rate is 0.12 m / day. The coverage area is 8 km × 8 km, with a resolution of 50 m × 50 m, totaling 25,600 grid points. Using a Bayesian update algorithm, this data is combined with historical geological datasets, including subsidence event records and geological parameters from the past five years, such as rock permeability of 0.0015 m / s and porosity of 0.25, to update the initial parameters in the numerical model. Next, a gradient descent optimization algorithm was employed with a learning rate of 0.01 and 100 iterations. The mean square error between the predicted values ​​and the actual monitoring data was calculated. For example, if the error at a certain grid point was 0.0004 m / day, which was below the set threshold of 0.0005 m / day, it indicated that the model parameters had effectively converged. Then, an adaptive Kalman filter algorithm was used to dynamically adjust the model parameters. Combined with real-time monitoring data (such as a groundwater level fluctuation standard deviation of 0.05 m at a certain point), predicted values ​​for the next three days were generated. For example, the predicted settlement rate for June 13, 2024, was 0.0032 m / day. Furthermore, using a spatial interpolation algorithm (Kriging), parameters for unmonitored areas were calculated based on the grid point data, resulting in a collapse probability exceeding 0.8 for approximately 1500 high-risk grid points. Finally, the optimized model parameters and prediction results were stored in NetCDF format and updated periodically to the geological disaster management platform via an automated script. This ensured that the model adapted to dynamic geological conditions, forming a closed-loop logical chain from spatiotemporal data integration to model optimization.

[0048] The above description is merely an example and illustration of the structure of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described, or use similar methods to replace them, as long as they do not deviate from the structure of the invention or exceed the scope defined in the claims, all of which should fall within the protection scope of the present invention.

Claims

1. A multi-media three-dimensional modeling and geotechnical mechanical response coupling method, characterized in that, The method includes: Data on groundwater level monitoring and geological exploration in karst areas are acquired. The monitoring data is then processed by a data preprocessing module to perform time-series standardization and spatial interpolation, resulting in a continuous spatiotemporal distribution of water level change sequences and geometric parameter information of the soil-rock interface. A numerical model of groundwater seepage field was established based on the water level change sequence. The fluid-structure interaction algorithm was used to calculate the pore water pressure distribution and seepage force field under different water level conditions, and the hydraulic gradient and seepage velocity distribution characteristics near the soil-rock interface were obtained. A stress field analysis model for soil and rock mass is constructed using the pore water pressure distribution data. If the water level rises, the top support stress effect is calculated; if the water level falls, the vacuum absorption stress effect is calculated, thus determining the stress concentration areas around the karst cave and on both sides of the karst conduit. The stress distribution data of the stress concentration area is used to calculate the soil deformation. The displacement field and strain field of the soil-rock interface are solved by the finite element method to obtain the spatial distribution and temporal evolution law of local deformation accumulation. A collapse risk assessment model is established based on the time-series evolution data of the accumulated deformation. If the local deformation rate exceeds a preset threshold, it is judged as a high-risk area. If the accumulated deformation reaches a critical value, the potential collapse location is determined. A spatial prediction model is constructed using the high-risk areas and potential collapse location information. A risk assessment algorithm is used to quantitatively calculate the collapse probability of different areas, and the spatial distribution prediction results of karst collapse are obtained. A dynamic early warning system is established based on the spatial distribution prediction results and real-time monitoring data. If the trend of the monitoring parameters changes is consistent with the output of the prediction model, the early warning mechanism is triggered to determine the time window and scope of impact of the collapse. The spatiotemporal information output by the aforementioned early warning mechanism is used to update the parameters of the numerical model, and the prediction accuracy is optimized and adjusted through a feedback correction mechanism to obtain a karst collapse prediction technology system adapted to complex geological conditions.

2. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process involves acquiring groundwater level monitoring data and geological exploration data in karst areas. A data preprocessing module performs time-series standardization and spatial interpolation on the monitoring data to obtain a continuous spatiotemporal distribution of water level change sequences and geometric parameter information of the soil-rock interface, including: The raw dataset is obtained by acquiring groundwater level monitoring data and geological exploration data from karst areas through a data acquisition system. A data preprocessing module was used to perform time-series standardization on the groundwater level monitoring data in the original dataset to generate water level change data with consistent time series. For water level change data, spatial interpolation methods are used to fill in missing areas and obtain a continuous spatiotemporal distribution of water level change sequence; Based on geological exploration data, the distribution characteristics of the soil-rock interface are extracted, and the geometric parameter information of the soil-rock interface is determined. If the water level change sequence shows abnormal fluctuations, then the geometric parameters of the soil-rock interface should be combined to determine whether there is a risk area for karst cave collapse. After obtaining the judgment results, based on the distribution of carbonate rocks and the characteristics of fracture development in the risk area, the spatial distribution pattern of the stratigraphic unconformity surface is analyzed to determine the potential collapse impact range; By analyzing the water level change trends and geological characteristics within the potential collapse area, a distribution layer of karst cave collapse risk is generated.

3. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process involves establishing a numerical model of the groundwater seepage field based on the water level change sequence, using a fluid-structure interaction algorithm to calculate the pore water pressure distribution and seepage force field under different water level conditions, and obtaining the hydraulic gradient and seepage velocity distribution characteristics near the soil-rock interface, including: Based on the water level change sequence and geological exploration data, a three-dimensional numerical model of the groundwater seepage field was constructed, and the pore water pressure distribution was calculated using the finite element method to obtain the seepage force field characteristics. By analyzing the characteristics of the seepage force field and combining the geometric parameters of the soil-rock interface, the hydraulic gradient and seepage velocity distribution are calculated to determine the seepage path within the overburden carbonate rock region. If the seepage velocity exceeds the preset threshold, the formation unconformity surface and fracture development characteristics are extracted to generate distribution data of high-risk seepage areas; Based on the distribution data of high-risk seepage areas, the spatial correlation between carbonate rock distribution and fracture network is analyzed to obtain the boundary information of the potential impact area of ​​cave collapse. By using the boundary information of the potential impact area and combining it with the hydraulic gradient distribution, a spatial distribution layer of collapse risk is generated using the Kriging interpolation method. Based on the spatial distribution layer of collapse risk, the superimposed areas of fracture development characteristics and stratigraphic unconformities are extracted to identify key areas for priority monitoring. By using monitoring data from key areas, the parameters of the seepage field model are updated to generate dynamically adjusted collapse risk distribution data.

4. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process involves constructing a stress field analysis model for the soil and rock mass using the pore water pressure distribution data. If the water level rises, the top-support stress effect is calculated; if the water level falls, the vacuum absorption stress effect is calculated. This determines the stress concentration areas around the karst cave and on both sides of the karst conduit, including: Based on the pore water pressure distribution data, a three-dimensional stress field model of rock and soil was constructed using the finite element method. The stress distribution characteristics under different water level changes were calculated, and the spatial distribution data of the stress concentration area was obtained. By using spatial distribution data of stress concentration areas, combined with the distribution characteristics of stratigraphic unconformities and fracture surfaces, spatial overlay analysis is employed to determine the overlapping areas of stress concentration areas and overlying carbonate rock distribution, thereby obtaining boundary information of high-risk areas. If the boundary information of the high-risk area overlaps with the distribution of karst conduits, then geostatistical analysis is used to calculate the stress gradient distribution within the stress concentration area to obtain the stress gradient anomaly area. Based on the stress gradient anomaly areas and combined with the geological features around the caves, the Kriging interpolation method is used to generate a spatial distribution layer of the stress gradient anomaly areas to identify potential collapse risk areas. By identifying potential collapse risk areas, the superimposed characteristics of fracture surface distribution and stratigraphic unconformity are extracted. Spatial clustering analysis is then used to obtain the priority ranking of high-risk areas and determine priority monitoring areas. Based on the stress distribution characteristics of the priority monitoring area, time series analysis is used to update the stress field model parameters and generate dynamically adjusted stress distribution data. By using dynamically adjusted stress distribution data, combined with the characteristics of top support stress and vacuum erosion stress, the finite difference method is used to calculate the stress change trend in high-risk areas, thus obtaining the dynamic distribution of potential collapse areas.

5. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process involves using stress distribution data from the stress concentration region to calculate soil deformation, and solving for the displacement and strain fields at the soil-rock interface using the finite element method. This yields the spatial distribution and temporal evolution of accumulated local deformation, including: Using the finite element method and based on stress distribution data, the displacement and strain fields at the soil-rock interface are calculated to obtain the spatial distribution characteristics and temporal evolution law of local deformation accumulation. Based on the spatial distribution characteristics of local deformation accumulation, combined with the planar distribution of karst cave collapse sections in the overburdened carbonate rock distribution area, spatial analysis methods are used to extract the superimposed characteristics of stratigraphic unconformity surfaces, fracture surfaces, and fracture zone development to identify potential high-risk areas for deformation. By combining the distribution data of the infill material with the potential high-risk deformation areas, the Kriging interpolation method is used to generate a spatial distribution layer of the deformation gradient within the high-risk areas, thereby obtaining the abnormal deformation gradient areas. If the abnormal deformation gradient area overlaps with the karst cave collapse area, then geostatistical analysis methods are used to calculate the deformation rate distribution within the abnormal area to obtain the abnormal deformation rate area. Based on the abnormal deformation rate areas, combined with the distribution characteristics of stratigraphic unconformities and fracture surfaces, spatial cluster analysis is used to determine the priority ranking of high-risk areas and obtain priority monitoring areas; By prioritizing the deformation rate distribution in the monitored area, time series analysis is used to update the displacement field and strain field model parameters, generating dynamically adjusted deformation distribution data. Based on the dynamically adjusted deformation distribution data, combined with the distribution of filling material and the development characteristics of the fracture zone, the finite difference method is used to calculate the deformation change trend of high-risk areas and determine the dynamic distribution of potential collapse areas.

6. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process of establishing a collapse risk assessment model based on the time-series evolution data of accumulated deformation, identifying a high-risk area if the local deformation rate exceeds a preset threshold, and determining the potential collapse location if the accumulated deformation reaches a critical value includes: By using time-series evolution data and employing time-series decomposition methods, the periodic and trend components of the local deformation rate are extracted to obtain the local deformation rate distribution. If the trend component of the local deformation rate exceeds a preset threshold, a spatial distribution layer of the high-risk area is generated by spatial interpolation to determine the range of the high-risk area. For high-risk areas, a threshold comparison method is used in conjunction with the cumulative deformation data. If the cumulative amount exceeds the critical value, a spatial positioning algorithm is used to determine the set of potential collapse locations. Based on the set of potential collapse locations and combined with data on carbonate rock distribution areas, spatial overlay analysis is used to extract the overlapping areas of stratigraphic unconformities and fracture surfaces to obtain geologically fragile areas. By analyzing geologically vulnerable areas and combining data on fracture zone characteristics, density analysis is used to calculate the fracture zone distribution density within the vulnerable areas and identify high-density fracture zone regions. For areas with high-density fractured zones, a classification algorithm is used to extract the distribution ratios of sand, clay, and gravel, based on data on the type of infill material, to obtain the distribution characteristics of the infill material. Based on the distribution characteristics of the filling material, a weighted overlay analysis was adopted to integrate data from areas with fragile geological structures and high-density fracture zones, generating a priority distribution of collapse risk and determining priority monitoring areas.

7. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The process involves constructing a spatial prediction model based on the high-risk areas and potential collapse locations, employing a risk assessment algorithm to quantitatively calculate the collapse probability in different areas, and obtaining spatial distribution prediction results for karst collapses, including: By using data on high-risk areas and potential collapse locations, and employing a kernel density estimation algorithm, the spatial probability distribution of karst collapse is calculated to obtain a collapse risk prediction layer. Based on the collapse risk prediction layer, the spatial range of the overburden carbonate rock distribution area is extracted. Combined with the distribution data of stratigraphic unconformities and fracture surfaces, spatial overlay analysis is used to determine the set of geologically vulnerable areas. If the set of geologically vulnerable areas overlaps with the distribution of fracture zones, spatial correlation analysis can be used to calculate the fracture zone density in the overlapping areas to obtain the high-density fracture zone areas. For areas with high-density fractured zones, data on the distribution of backfill material was obtained. A classification algorithm was used to extract the distribution ratios of sand, clay, and gravel, thus obtaining the characteristic distribution of the backfill material. Based on the distribution characteristics of the infill material and combined with the set of geologically vulnerable areas, a weighted overlay analysis was used to generate a collapse sensitivity distribution layer. By using a collapse sensitivity distribution layer and employing a spatial clustering algorithm, the set of monitoring areas with the highest priority is determined. For the highest priority monitoring area set, combined with historical collapse event data, time series analysis is used to predict the probability distribution of future collapses, thus obtaining a dynamic risk distribution.

8. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The dynamic early warning system established based on the spatial distribution prediction results and real-time monitoring data triggers an early warning mechanism if the trend of monitoring parameter changes is consistent with the output of the prediction model, determining the time window and impact range of the collapse, including: By using spatial distribution data and prediction results, stratigraphic distribution information within the overburden carbonate rock area is obtained. Combined with real-time monitoring data, abnormal areas of parameter changes are identified, and a preliminary risk distribution range is obtained. Based on the preliminary risk distribution range, relevant data on fracture fracturing are obtained for abnormal areas with parameter changes. Spatial overlay analysis is used to determine the degree of overlap between fracture fracturing and stratigraphic distribution, and to identify the set of high-risk areas. If the parameter changes within the high-risk area set exceed the preset threshold, the change trend within the time window is extracted by real-time monitoring data to obtain the distribution of potential collapse periods. Based on the distribution of potential collapse periods and combined with the data on the scope of impact, the distribution characteristics of the filling material are obtained. A classification algorithm is used to determine the weight of the impact of the filling material on the collapse, and a risk-weighted distribution is obtained. By using risk-weighted distribution, real-time monitoring update data is obtained to meet the needs of dynamic early warning. If the update data is consistent with the trend of risk-weighted distribution, the early warning mechanism is triggered to determine the priority early warning area. Based on the priority warning areas, combined with stratigraphic distribution and fracture data, a spatial clustering algorithm is used to extract key monitoring points within the affected area, thus obtaining a set of key monitoring objects; By focusing on a set of key monitoring targets, a continuous real-time data stream is obtained. If the continuous data stream shows that parameter changes are consistently abnormal, the predicted value of the time window is updated to determine the final warning time range.

9. The multi-media three-dimensional modeling and geotechnical mechanical response coupling method according to claim 1, characterized in that, The numerical model parameters are updated using the spatiotemporal information output by the aforementioned early warning mechanism, and the prediction accuracy is optimized and adjusted through a feedback correction mechanism to obtain a karst collapse prediction technology system adapted to complex geological conditions, including: The spatiotemporal information obtained through the early warning mechanism is used to dynamically adjust the parameters of the numerical model. Information processing is performed on the data characteristics under complex geological environments to obtain an updated set of model parameters. Based on the updated model parameter set, a feedback correction mechanism is used to calibrate the prediction accuracy. Combining planar distribution data and the regional characteristics of carbonate rocks, the range of prediction deviation after calibration is determined. If the calibrated prediction deviation exceeds the preset threshold, the distribution characteristics of potential karst collapse areas are extracted through the fusion processing of spatiotemporal information and geological environment data to obtain the high-risk area classification results. Based on the high-risk area delineation results, data on fracture distribution and fracture zone areas were obtained. Spatial analysis tools were used to compare the degree of overlap between fracture distribution and planar distribution to determine the set of potential collapse areas. If the fracture distribution characteristics within the potential collapse area set are highly correlated with the data of the fracture zone area, then the distribution information of the filling material can be used for auxiliary analysis, combined with the geological environment characteristics, to determine the scope of key areas of concern. Based on the key areas of concern, the dynamic changes of karst collapse are monitored using real-time data streams. Combined with the output of the feedback correction mechanism, an updated risk distribution layer is obtained. By using the updated risk distribution layer and combining it with the spatiotemporal information output of the early warning mechanism, a logistic regression algorithm is used to classify the risk distribution layer and determine the final set of high-risk warning areas.

Images