A method for assessing ground subsidence risk through multi-source information fusion

Abstract

This application relates to a multi-source information fusion method for ground subsidence risk assessment, belonging to the field of geological monitoring technology. It involves collecting static geological environment data and dynamic environmental disturbance data, performing unified spatiotemporal benchmark conversion and standardization; constructing a multi-dimensional ground subsidence risk assessment factor database based on standard grid units; extracting spatial and temporal features to construct a risk assessment sample set; employing a CNN-LSTM hybrid neural network model to simultaneously extract spatial neighborhood features and temporal evolution features; using a generative adversarial network for sample balancing and hyperparameter optimization; and outputting the risk assessment results for visualization. This invention significantly improves the fusion capability of multi-source heterogeneous data, employs a CNN-LSTM hybrid model to simultaneously extract spatial and temporal features, significantly improves the prediction accuracy of the risk assessment model, effectively solves the problem of scarce high-risk samples, and significantly enhances the accuracy and reliability of ground subsidence risk assessment.

Description

Technical Field

[0001] This application relates to a ground subsidence risk assessment method based on multi-source information fusion, belonging to the field of geological monitoring technology. Background Technology

[0002] With the continuous advancement of urbanization and the large-scale development of urban underground space, ground subsidence has become a significant type of geological hazard affecting urban safety. Urban ground subsidence is typically triggered by a variety of factors, including changes in groundwater levels, leakage from underground pipelines, karst development, subsidence of mining subsidence areas, and disturbance during foundation pit construction. Its occurrence is characterized by its insidious, sudden, and highly random nature, posing a serious threat to the safe operation of urban roads, rail transit, underground pipelines, and surrounding buildings and structures. Therefore, establishing a scientific and effective method for assessing urban ground subsidence risk is of great significance for urban geological hazard prevention and control, underground space development safety, and emergency decision-making.

[0003] The methods for assessing urban ground subsidence risk have evolved from qualitative judgment to quantitative evaluation. Early assessments primarily relied on geological survey data, historical disaster records, and expert experience, using qualitative analysis of regional geological conditions to determine risk levels. Subsequently, risk zoning methods based on evaluation index systems and GIS spatial overlay analysis gradually developed. These methods typically first collect influencing factors such as topography, stratigraphy, groundwater conditions, road networks, historical subsidence points, land use, and engineering activities. Then, each factor is graded and assigned a value. Weights are then determined using the analytic hierarchy process (AHP), entropy weighting, expert scoring, or a combination of weighting methods. Finally, the factor layers are overlaid and analyzed in a GIS platform to obtain a risk zoning map. This type of method has a clear workflow, is easy to implement, and can complete regional-scale risk assessments based on existing geological and historical disaster data.

[0004] However, existing technologies typically separate static geological condition evaluation from dynamic monitoring and analysis, making it difficult to achieve unified fusion and collaborative modeling of multi-source heterogeneous data. Traditional GIS evaluation methods focus on the spatial overlay of static factors, failing to effectively utilize time-series information from dynamic monitoring data such as InSAR deformation, rainfall, and groundwater levels. While machine learning methods can learn the nonlinear relationships between multiple factors, their input features are mostly limited to static factors, lacking the ability to express temporal evolution processes. Although time-series prediction models can capture the development trend of subsidence or deformation, they do not incorporate static background conditions such as geological structures, borehole parameters, and hydrogeology into a unified framework. Furthermore, urban ground subsidence is a low-frequency disaster event, with a limited number of historical subsidence samples and a low proportion of high-risk samples, making model training susceptible to imbalanced samples. These issues make it difficult for existing technologies to simultaneously complete multi-source data fusion, spatial feature extraction, temporal feature analysis, and dynamic risk expression within the same framework, requiring improvements in the accuracy, stability, and engineering applicability of the evaluation results. Summary of the Invention

[0005] The purpose of this invention is to provide an intelligent evaluation method for urban ground subsidence, which overcomes the shortcomings of existing ground subsidence risk assessment technologies, such as insufficient fusion of multi-source heterogeneous data, weak collaborative analysis capability of static geological conditions and dynamic triggering factors, difficulty in unified modeling of spatial neighborhood characteristics and temporal evolution characteristics, small number of high-risk samples and unbalanced categories, and insufficient dynamic expression capability of evaluation results.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a ground subsidence risk assessment method based on multi-source information fusion, comprising the following steps:

[0007] Step 1: Collect multi-source data required for urban ground subsidence risk assessment. The multi-source data includes static geological environment data and dynamic environmental disturbance data.

[0008] Step 2: Perform unified spatiotemporal benchmark transformation and standardization on the multi-source data to unify it to a unified spatial grid;

[0009] Step 3: Construct a multi-dimensional risk assessment factor database for ground subsidence based on the unified spatial grid;

[0010] Step 4: Extract spatial features from static factors and extract temporal features from dynamic factors from the risk assessment factor database to construct a risk assessment sample set;

[0011] Step 5: Construct a CNN-LSTM ground subsidence risk assessment model based on the risk assessment sample set. The model includes a spatial feature extraction branch, a temporal feature extraction branch, a feature fusion layer, and a risk output layer.

[0012] Step 6: The CNN-LSTM ground collapse risk assessment model is optimized by sample balancing and model hyperparameter optimization, and the risk assessment results are output and visualized.

[0013] Furthermore, the static geological environment data mentioned in step 1 includes:

[0014] Topographic and geomorphological data, stratigraphic and lithological data, geological structure data, borehole data, hydrogeological data, surface cover data, underground space development data, engineering survey data, road and pipeline network data, and historical ground subsidence disaster data;

[0015] The dynamic environmental disturbance data includes InSAR surface deformation monitoring data, rainfall data, groundwater level monitoring data, engineering construction disturbance data, surface subsidence rate data, and cumulative subsidence data.

[0016] Furthermore, in step 2, the multi-source data undergoes unified coordinate transformation, projection processing, pruning, resampling, outlier removal, missing value repair, and field standardization; thus unifying all types of data to the same spatial reference system and a unified time scale.

[0017] Furthermore, the database mentioned in step 3 includes:

[0018] Static factor layer, dynamic factor layer, sample label layer, and model output layer;

[0019] The static factor layer uses a standard grid as the basic unit to store topographic and geomorphological factors, stratigraphic and lithological factors, hydrogeological factors, geological structural factors, borehole interpolation factors, land cover factors, historical disaster proximity factors, and engineering disturbance background factors.

[0020] The dynamic factor layer uses grid number and time number as a joint index to store cumulative surface subsidence, subsidence rate over time, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change, and rainfall lag characteristics.

[0021] Furthermore, in step 4, when extracting spatial features from static factors, the values ​​of similar factors within the surrounding neighborhood of the target grid are extracted as spatial neighborhood features.

[0022] When extracting time-series features of dynamic factors, a time window is constructed with each grid cell as the object to extract cumulative settlement, average settlement rate, maximum settlement rate, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change amplitude, and rainfall lag term.

[0023] Furthermore, the spatial feature extraction branch in step 5 employs a convolutional neural network to extract the spatial structure features of the target grid and its neighborhood.

[0024] The time feature extraction branch uses a long short-term memory network to extract the time evolution features of dynamic factors;

[0025] The feature fusion layer fuses spatial and temporal features; the risk output layer outputs the probability or level of ground subsidence risk.

[0026] Furthermore, the sample balancing process in step 6 employs at least one of the following methods: undersampling, oversampling, class weight adjustment, and generative adversarial network enhancement.

[0027] The hyperparameter optimization employs an automated search method to optimize the spatial neighborhood window size, the number of convolutional layers, the number of convolutional kernels, the number of LSTM hidden units, the learning rate, and the batch size.

[0028] Furthermore, in step 6, areas with a risk probability less than the first preset value are classified as low-risk areas, areas with a risk probability greater than the first preset value and less than the second preset value are classified as medium-risk areas, areas with a risk probability greater than the second preset value and less than the third preset value are classified as relatively high-risk areas, and areas with a risk probability greater than the third preset value are classified as high-risk areas.

[0029] Furthermore, the visualization described in step 6 includes:

[0030] Risk probability grid map, risk level zoning map, list of key risk grids, risk factor contribution results, time series risk change curves and model uncertainty results;

[0031] The risk assessment results will be integrated into the geographic information system platform and the 3D visualization platform for comprehensive display.

[0032] Furthermore, the method also includes performing uncertainty quantification analysis on key monitoring points and evaluating the stability of the points by calculating the standard deviation of time-series fluctuations;

[0033] When the standard deviation of the time series fluctuation is greater than a preset threshold, the region is marked as a high uncertainty region and a prompt is made to prioritize manual on-site verification.

[0034] Compared with the prior art, the present invention has the following beneficial effects:

[0035] First, by unifying the coordinate system, spatial grid, time scale, and field standards, this invention incorporates data such as geological structure, stratigraphy, borehole data, hydrogeology, historical subsidence points, InSAR deformation, rainfall, and groundwater level into the same database, achieving unified management of static geological background and dynamic triggering factors, and significantly improving the fusion capability of multi-source heterogeneous data.

[0036] Second, this invention transforms discrete borehole parameters into a continuous geological parameter field by using methods such as inverse distance weighting, Kriging interpolation, and Gaussian process regression, thereby improving the continuity and integrity of the expression of underground geological conditions and significantly improving the accuracy of spatial reconstruction of geological parameters.

[0037] Third, this invention uses a CNN-LSTM hybrid model to simultaneously extract spatial neighborhood features and time series features. The CNN branch learns the spatial distribution patterns of geological structures, strata lithology, borehole interpolation parameters, and historical disaster density, while the LSTM branch learns the dynamic temporal patterns of InSAR deformation, rainfall, and groundwater level changes. The fusion of the two significantly improves the prediction accuracy of the ground subsidence risk assessment model.

[0038] Fourth, this invention improves the model's ability to identify high-risk samples by using generative adversarial networks to enhance generative samples, adjust class weights, and optimize thresholds, effectively solving the problems of a small number of urban ground subsidence event samples and a low proportion of high-risk samples.

[0039] Fifth, this invention further converts the continuous risk score output by the model into low risk, medium risk, relatively high risk and high risk levels, and determines the risk discrimination criteria through threshold optimization. At the same time, it introduces risk confidence, uncertainty analysis and spatial consistency verification, which significantly improves the operability and reliability of the risk zoning results.

[0040] Sixth, this invention outputs risk results through two-dimensional GIS and three-dimensional visualization, enabling the evaluation results to be directly used for risk zoning, early warning analysis, engineering management and decision support, significantly improving the efficiency of result expression and engineering application. Attached Figure Description

[0041] Figure 1 This is a flowchart illustrating the process of a ground subsidence risk assessment method based on multi-source information fusion, as claimed in an embodiment of the present invention. Detailed Implementation

[0042] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of this application.

[0043] The terms "first," "second," and "third" in this application are for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first," "second," or "third" may explicitly or implicitly include at least one of those features. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified. All directional indications in the embodiments of this application, such as up, down, left, right, front, back, etc., are only used to explain the relative positional relationships and movements between components in a specific orientation as shown in the accompanying drawings. If the specific orientation changes, the directional indications will change accordingly. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices.

[0044] References to embodiments herein mean that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a mutually exclusive, independent, or alternative embodiment. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0045] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0046] According to the first embodiment of the present invention, referring to Figure 1 This invention claims protection for a ground subsidence risk assessment method based on multi-source information fusion, comprising the following steps:

[0047] Step 1: Collect multi-source data required for urban ground subsidence risk assessment. The multi-source data includes static geological environment data and dynamic environmental disturbance data.

[0048] Step 2: Perform unified spatiotemporal benchmark transformation and standardization on the multi-source data to unify it to a unified spatial grid;

[0049] Step 3: Construct a multi-dimensional risk assessment factor database for ground subsidence based on the unified spatial grid;

[0050] Step 4: Extract spatial features from static factors and extract temporal features from dynamic factors from the risk assessment factor database to construct a risk assessment sample set;

[0051] Step 5: Construct a CNN-LSTM ground subsidence risk assessment model based on the risk assessment sample set. The model includes a spatial feature extraction branch, a temporal feature extraction branch, a feature fusion layer, and a risk output layer.

[0052] Step 6: The CNN-LSTM ground collapse risk assessment model is optimized by sample balancing and model hyperparameter optimization, and the risk assessment results are output and visualized.

[0053] Furthermore, the static geological environment data mentioned in step 1 includes:

[0054] Topographic and geomorphological data, stratigraphic and lithological data, geological structure data, borehole data, hydrogeological data, surface cover data, underground space development data, engineering survey data, road and pipeline network data, and historical ground subsidence disaster data;

[0055] The dynamic environmental disturbance data includes InSAR surface deformation monitoring data, rainfall data, groundwater level monitoring data, engineering construction disturbance data, surface subsidence rate data, and cumulative subsidence data.

[0056] In this embodiment, step 1 involves collecting multi-source data required for urban ground subsidence risk assessment, specifically including two main categories: static geological environment data and dynamic environmental disturbance data. The sources of static geological environment data include topographic and geomorphological data, stratigraphic and lithological data, geological structure data, borehole data, hydrogeological data, land cover data, underground space development data, engineering survey data, road and pipeline network data, and historical ground subsidence disaster point data. In actual implementation, static geological environment data includes at least six categories of geological safety influencing factors, covering topographic and geomorphological factors, hydrogeological factors, land cover factors, geological structure factors, borehole parameter factors, and historical disaster background factors. This data can be obtained from geological survey departments, urban surveying and mapping departments, underground pipeline management departments, geological disaster monitoring agencies, and historical disaster archives. The static data acquisition formats include vector shapefile format, raster GeoTIFF format, CAD drawing format, and text table format. The coordinate system of the data source unit may be the Beijing 54 coordinate system, the Xi'an 80 coordinate system, or the National 2000 coordinate system.

[0057] The sources of dynamic environmental disturbance data include InSAR surface deformation monitoring data, rainfall data, groundwater level monitoring data, engineering construction disturbance data, surface subsidence rate data, and cumulative subsidence data. InSAR surface deformation monitoring data is acquired through synthetic aperture radar satellites, including Sentinel-1, ALOS-2, and Radarsat-2, with a spatial resolution of 5 to 25 meters, a sampling period of 12 days, and a single image covering tens of thousands of square kilometers. Rainfall data is acquired through automatic observations at meteorological stations, including daily rainfall, hourly rainfall, and maximum hourly rainfall intensity. Data sources include national meteorological station networks, urban meteorological observation network networks, and regional automatic weather station networks. Groundwater level monitoring data is obtained through monitoring wells, located in key subsidence areas, along subway lines, underground engineering construction areas, and areas surrounding historical subsidence sites. Engineering construction disturbance data includes excavation depth, tunnel boring machine (TBM) parameters, and pipe jacking operation records, which can be obtained from construction units or engineering supervision units.

[0058] Furthermore, in step 2, the multi-source data undergoes unified coordinate transformation, projection processing, pruning, resampling, outlier removal, missing value repair, and field standardization; thus unifying all types of data to the same spatial reference system and a unified time scale.

[0059] In this embodiment, step 2 involves performing a unified spatiotemporal benchmark transformation and standardization on the multi-source data collected in step 1. The unified coordinate transformation process unifies data from different sources to the WGS-84 coordinate system and the UTM projection system. The coordinate transformation method uses a seven-parameter transformation model or a four-parameter transformation model. The transformation parameters are obtained through control points, and the control points are at least level three or lower level points. The projection processing process converts coordinates in the geographic coordinate system to planar coordinates in the projected coordinate system. The selection of the projection zone is determined based on the longitude range of the study area, typically choosing the projection zone whose central meridian longitude is closest to the center longitude of the study area. The clipping process clips the data according to the administrative boundaries or custom boundaries of the study area. The boundary range is defined using polygon files in shapefile format, and the clipping tool uses the GDAL library or the Clip tool on the ArcGIS platform. The resampling process uniformly resamples all types of raster data to a standard 30m x 30m grid. The resampling method is determined based on the data type; for continuous variables, bilinear interpolation or cubic convolution interpolation is used, and for categorical variables, the nearest neighbor method is used. The resampled data is stored in GeoTIFF format, with numerical precision retained to six decimal places.

[0060] The outlier removal process identifies and processes clearly unreasonable values ​​in various data types. Outliers in InSAR deformation data are characterized by significant scatterer loss, atmospheric delay phase remnants, or systematic biases caused by orbital errors. Removal methods employ temporal consistency checks and spatial neighborhood comparisons, with a threshold set at three times the standard deviation of the deformation rate. Outliers in precipitation data are characterized by negative rainfall, extremely high rainfall exceeding historical extremes, or abrupt changes after prolonged periods without rainfall. Removal methods employ comparisons with neighboring stations and historical data for the same period. Outliers in groundwater level data are characterized by abrupt increases or decreases, exceeding aquifer limits, or changes that clearly violate hydrodynamic principles. Removal methods employ water level fluctuation threshold control and correlation checks with adjacent monitoring wells.

[0061] The missing value repair process involves interpolating to fill in gaps in the data. Incomplete records in borehole data are repaired using inverse distance weighted interpolation or Kriging interpolation, with 8 to 20 neighboring boreholes involved in the interpolation. The interpolation search radius is determined based on borehole density, typically not exceeding 5 kilometers. Missing sections in rainfall data are filled using linear or spline interpolation, with the filling period not exceeding three consecutive observation cycles; missing sections exceeding this range are marked as invalid data. Missing sections in groundwater level data are filled using time series fitting methods, including polynomial fitting or Fourier series fitting, with the fitting order determined based on the periodic characteristics of water level changes. Isolated noise points in InSAR results are removed using median filtering or spatial neighborhood filtering, with the filtering window size set to 3x3 to 5x5 pixels.

[0062] The field standardization process organizes data from different sources according to unified format requirements. Attribute field names uniformly adopt English abbreviations, and the naming rules follow the requirements of the "Geographic Information Element Naming Specification". Data types are uniformly set to numeric, text, or date. Numeric fields use international standard units, text fields use a unified encoding table for value assignment, and date fields uniformly adopt the ISO 8601 standard format.

[0063] Furthermore, the database mentioned in step 3 includes:

[0064] Static factor layer, dynamic factor layer, sample label layer, and model output layer;

[0065] The static factor layer uses a standard grid as the basic unit to store topographic and geomorphological factors, stratigraphic and lithological factors, hydrogeological factors, geological structural factors, borehole interpolation factors, land cover factors, historical disaster proximity factors, and engineering disturbance background factors.

[0066] The dynamic factor layer uses grid number and time number as a joint index to store cumulative surface subsidence, subsidence rate over time, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change, and rainfall lag characteristics.

[0067] In this embodiment, step 3 involves constructing a multi-dimensional ground subsidence risk assessment factor database based on a unified spatial grid. The database adopts a hierarchical storage structure, including a static factor layer, a dynamic factor layer, a sample label layer, and a model output layer. The static factor layer uses a standard grid as the basic unit, with each grid unit assigned a unique CellID. The CellID is encoded using a combination of row and column numbers, starting from 1 and sequentially. The total number of rows and columns is calculated based on the study area and grid size. The static factor layer stores factor types including topographic factors, stratigraphic lithology factors, hydrogeological factors, geological structure factors, borehole interpolation factors, land cover factors, historical disaster proximity factors, and engineering disturbance background factors. A total of 14 categories of specific parameters are set for each factor.

[0068] Topographic factors include elevation, slope, aspect, and topographic relief. Elevation data is derived from a digital elevation model (DEM) with a resolution of 30 meters; slope is calculated from DEM data using the maximum gradient method; aspect is calculated from DEM data using the differential method; topographic relief is obtained through local window statistics, with a window size of 11 x 11 pixels. Stratigraphic lithology factors are stored in stratigraphic code form, using the stratigraphic division scheme specified in the "Regional Geological Survey Specification". Hydrogeological factors include groundwater depth, aquifer thickness, and permeability coefficient. Geological structural factors include fault zone distance and active fault markers, with distances calculated using Euclidean distance or network distance. Borehole interpolation factors are obtained through interpolation processing in step 4, including water level depth interpolation, weak layer thickness interpolation, and sand layer thickness interpolation. Land cover factors are stored using land use type classification codes, using the national land resource classification standard; historical disaster proximity factors are calculated using kernel density estimation.

[0069] The dynamic factor layer organizes data using a combined index of CellID and TimeID. TimeID is encoded as year-month-day. The dynamic factor layer stores factor types including cumulative surface subsidence, time-period subsidence rate, cumulative rainfall, maximum daily rainfall, API (Advanced Precipitation Index), groundwater level change, and rainfall lag characteristics from 1 to 7 days. A total of 8 parameter categories are set for each factor. Cumulative surface subsidence is calculated as the cumulative InSAR deformation value within the time window. Time-period subsidence rate is calculated by dividing the deformation value within the time window by the time window length. Cumulative rainfall is calculated as the cumulative daily rainfall value within the time window. The API is calculated by summing the attenuated daily rainfall from the previous period, with an attenuation coefficient set to 0.9. Rainfall lag characteristics include 1-day, 3-day, 5-day, and 7-day delayed rainfall, reflecting the triggering effect of rainfall with different lag periods on ground subsidence.

[0070] The database employs a hierarchical storage approach. Raster data is stored in GeoTIFF format, with single-band storage and compression methods including LZW or DEFLATE. Vector data is stored in Shapefile format, encompassing point, line, and polygon features, with spatial index files generated synchronously. Tabular data is stored in CSV or Parquet format, with commas as field separators and UTF-8 character encoding. Relationships between different data layers are established using three key fields: CellID, TimeID, and PointID. Relational queries utilize database table joins or the `Merge` function from the pandas library in Python.

[0071] Furthermore, in step 4, when extracting spatial features from static factors, the values ​​of similar factors within the surrounding neighborhood of the target grid are extracted as spatial neighborhood features.

[0072] When extracting time-series features of dynamic factors, a time window is constructed with each grid cell as the object to extract cumulative settlement, average settlement rate, maximum settlement rate, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change amplitude, and rainfall lag term.

[0073] In this embodiment, step 4 involves spatial feature extraction of static factors, transforming various factors from their original data format into a standard grid feature matrix. The dimension of the feature matrix is ​​the total number of grids multiplied by the feature dimension. The total number of grids is calculated based on the study area and grid size, while the feature dimension equals the number of static factors plus the spatial neighborhood feature dimension. Spatial neighborhood feature extraction uses the target grid as the center, extracting similar factor values ​​within a 5x5 or 7x7 grid radius as convolution input. For topographic factors, the mean, standard deviation, extreme values, and texture features within the neighborhood are extracted; for stratigraphic and lithological factors, the dominant lithology code and lithology conversion frequency within the neighborhood are extracted; and for hydrogeological factors, the water level gradient and hydraulic connection strength within the neighborhood are extracted.

[0074] For borehole data processing, extracted parameters include water level depth, borehole depth, inflow rate, flow rate, drawdown, soil layer thickness, weak layer thickness, sand layer thickness, and clay layer thickness. Spatial interpolation methods employed include the inverse distance weighting method, ordinary kriging, radial basis function method, or Gaussian process regression. For the inverse distance weighting method, the power exponent was set to 2, the search radius was set to a maximum influence distance of 5 km, and the number of neighboring points was set to 12. For the ordinary kriging method, a spherical or exponential model was selected for the variogram, the nugget value was set to 0, and the sill and range values ​​were determined based on the experimental variogram fitting. For the radial basis function method, multiple quadratic splines or thin-plate splines were selected as the basis functions. For the Gaussian process regression method, a radial basis function kernel or a Matrn kernel was selected as the kernel function, and the hyperparameters were determined through maximum likelihood estimation. The number of neighboring boreholes involved in the calculation for each interpolation point ranged from 8 to 20, and the interpolation results were output in raster form with a resolution consistent with the standard grid.

[0075] Temporal features were extracted from dynamic factors. A time window was constructed for each CellID or InSAR monitoring point, with a time step length of 6 to 12 time steps. Extracted feature types included cumulative settlement, average settlement rate, maximum settlement rate, cumulative rainfall, maximum daily rainfall, antecedent rainfall index (API), groundwater level variation, and rainfall lag term. Cumulative settlement was calculated as the sum of InSAR deformation values ​​at all time points within the time window. The average settlement rate was calculated by dividing the cumulative settlement by the time window length. The results showed that the maximum settlement rate reflected the maximum deformation rate at a single point within the time window. Cumulative rainfall was calculated as the sum of rainfall values ​​for all observed days within the time window. Maximum daily rainfall was calculated as the maximum rainfall value for all observed days within the time window. The antecedent rainfall index (API) was calculated by weighting and summing antecedent rainfall using an exponential decay function. The groundwater level variation was calculated as the difference between the highest and lowest groundwater levels within the time window. The rainfall lag term is constructed by correlating the rainfall at the current time point with the rainfall at the previous 1, 3, 5, and 7 time points, respectively.

[0076] After feature extraction, continuous variables are standardized. For normally distributed continuous variables, Z-score standardization is used, with the formula being the variable value minus the mean and then divided by the standard deviation. For non-normally distributed continuous variables, Min-Max normalization is used, with the formula being the variable value minus the minimum and then divided by the difference between the maximum and minimum values. Standardization parameters are calculated based on the training set, and the same standardization parameters are used for both the test and prediction sets. Categorical variables are encoded, either by one-hot encoding to convert them into binary vectors or by label encoding to convert them into integer indices.

[0077] The labeling of the risk assessment sample set is based on the spatial overlay relationship between historical ground subsidence disaster sites and standard grids. High-risk samples are identified by labeling the grid containing the historical subsidence disaster site and the surrounding buffer zone as high-risk, with the buffer radius set to 30 to 100 meters. Low-risk samples are identified by randomly selecting grid cells in areas far from historical subsidence sites, with low deformation rates and no obvious adverse geological conditions, and labeling them as low-risk. The sample set is partitioned using spatial grouping, dividing the study area into several spatial sub-regions based on geographical location. Samples within each sub-region are proportionally allocated to the training, validation, and test sets, with ratios of 6:2:2, 7:2:1, or 8:1:1. Spatial grouping effectively reduces data leakage caused by spatial autocorrelation and improves the model's generalization ability in unknown areas.

[0078] Furthermore, the spatial feature extraction branch in step 5 employs a convolutional neural network to extract the spatial structure features of the target grid and its neighborhood.

[0079] The time feature extraction branch uses a long short-term memory network to extract the time evolution features of dynamic factors;

[0080] The feature fusion layer fuses spatial and temporal features; the risk output layer outputs the probability or level of ground subsidence risk.

[0081] In this embodiment, step 5 involves constructing a CNN-LSTM ground subsidence risk assessment model. The overall architecture comprises four components: a spatial feature extraction branch, a temporal feature extraction branch, a feature fusion layer, and a risk output layer. The spatial feature extraction branch is implemented using a convolutional neural network. The input data is a multi-channel spatial feature tensor, with the tensor's dimension being the batch size multiplied by the number of channels, the neighborhood window height, and the neighborhood window width. The neighborhood window size is set to 5x5 or 7x7 grid cells, and the number of channels equals the number of static factors. The number of convolutional layers is set to 2 to 4, with each convolutional layer followed by a batch normalization layer and an activation function layer. The convolutional kernel size is set to 3x3 or 5x5, and the number of convolutional kernels is set to 32 to 128. The activation function is either ReLU or LeakyReLU. The pooling layer uses max pooling or average pooling, with a pooling window size of 2x2 and a stride of 2.

[0082] The temporal feature extraction branch employs a Long Short-Term Memory (LSTM) network. The input data is a time-series feature matrix, with dimensions equal to the batch size multiplied by the time step length and the number of dynamic factors. The time step length is set to 6 to 12 time steps, and the number of dynamic factors is 8 classes. The number of LSTM hidden units is set to 32 to 128, and the number of LSTM layers is set to 1 to 2. The Dropout ratio is set to 0.1 to 0.3 to prevent overfitting. A bidirectional LSTM can be optionally enabled as needed; when enabled, it can simultaneously learn forward and backward temporal dependencies.

[0083] The feature fusion layer fuses the spatial feature vector output from the spatial feature extraction branch and the temporal feature vector output from the temporal feature extraction branch. Fusion methods can include concatenation fusion, weighted fusion, or attention fusion. Concatenation fusion concatenates the two feature vectors dimensionally, resulting in a vector whose dimension is the sum of the spatial and temporal feature dimensions. Weighted fusion multiplies each feature vector by a learnable weight coefficient, which is then summed. Attention fusion introduces an attention mechanism, calculating attention weights to weight the spatial and temporal features. These attention weights can be calculated using scaled dot product attention or multi-head attention mechanisms. Following the fusion layer is a fully connected layer with 16 to 128 neurons. The output of the fully connected layer serves as the input to the risk output layer.

[0084] The risk output layer outputs risk probability or risk level based on task requirements. When outputting risk probability, a Sigmoid activation function is used, mapping the output value to the interval between 0 and 1, representing the probability of ground subsidence. When outputting multi-level risk levels, a Softmax activation function is used, outputting the probability distributions for four categories: low risk, medium risk, relatively high risk, and high risk. The choice of model loss function is determined based on the task type; binary cross-entropy loss is used for binary classification tasks, and multi-class cross-entropy loss is used for multi-class classification tasks. For cases where the number of high-risk samples is small, a weighted cross-entropy loss function or Focal Loss loss function is used, increasing the loss weight of the high-risk category to improve the model's ability to identify high-risk areas.

[0085] The overall training process of the model includes four stages: forward propagation, loss calculation, backpropagation, and parameter update. In the forward propagation process, the input data is sequentially passed through a spatial feature extraction branch, a temporal feature extraction branch, a feature fusion layer, and a risk output layer to obtain the predicted output. The loss calculation process compares the predicted output with the true labels and calculates the loss function value. The backpropagation process calculates the gradient of the parameters of each layer based on the loss function value. The parameter update process uses an optimization algorithm to update the model parameters; the optimization algorithm can be Adam, SGD, or AdamW. The learning rate is set to 0.0001 to 0.001, and the batch size is set to 16 to 64.

[0086] Furthermore, the sample balancing process in step 6 employs at least one of the following methods: undersampling, oversampling, class weight adjustment, and generative adversarial network enhancement.

[0087] The hyperparameter optimization employs an automated search method to optimize the spatial neighborhood window size, the number of convolutional layers, the number of convolutional kernels, the number of LSTM hidden units, the learning rate, and the batch size.

[0088] Furthermore, in step 6, areas with a risk probability less than the first preset value are classified as low-risk areas, areas with a risk probability greater than the first preset value and less than the second preset value are classified as medium-risk areas, areas with a risk probability greater than the second preset value and less than the third preset value are classified as relatively high-risk areas, and areas with a risk probability greater than the third preset value are classified as high-risk areas.

[0089] In this embodiment, step 6 addresses the issues of a small number of historical ground subsidence samples and a low proportion of high-risk samples by employing sample balancing and generative data augmentation methods to improve the model's learning ability for high-risk samples. Sample balancing methods include undersampling, oversampling, class weight adjustment, and Generative Adversarial Network (GAN) enhancement. Undersampling randomly removes some samples from the majority class to bring the number of majority and minority class samples closer to balance. Oversampling repeatedly samples minority class samples to increase their number. Class weight adjustment sets different weight coefficients for different classes in the loss function, with the weight coefficient for high-risk classes set to 2 to 5 times that of low-risk classes.

[0090] The Generative Adversarial Network (GAN) enhancement method generates synthetic high-risk samples by training a generator and a discriminator. The generator's input consists of a random noise vector and a risk class conditional vector, and its output is the feature vector of the synthetic high-risk sample. The random noise vector has a dimension of 32 to 128, and the conditional vector is a one-hot encoded risk class label. The network structure of the generator and discriminator can be a fully connected network or a convolutional network. The fully connected network has 2 to 5 layers with 64 to 512 neurons per layer. The convolutional network has 3 to 6 layers with 32 to 128 kernels per layer. The GAN is trained for 500 to 2000 epochs, with a batch size of 32 to 128 and a learning rate of 0.0001 to 0.001. During training, the generator and discriminator are updated alternately. The generator's optimization objective is to maximize the discriminator's misclassification probability of synthetic samples, while the discriminator's optimization objective is to maximize the correct classification probability of real and synthetic samples.

[0091] The sample selection process is evaluated based on four dimensions. Feature distribution similarity assessment uses principal component analysis or t-SNE to project synthetic and real samples into a low-dimensional space and calculate their distribution similarity. Temporal trend consistency assessment checks whether the time-series patterns of the synthetic samples conform to the evolutionary laws of ground subsidence. Spatial location rationality assessment checks whether the geographic coordinates of the synthetic samples are located within a reasonable study area. Geological mechanism rationality assessment checks whether the attribute values ​​of the synthetic samples conform to geological laws; for example, the thickness of weak layers should not exceed a reasonable range, and the groundwater level depth should not be negative. The synthetic samples selected through this four-dimensional evaluation, together with the real samples, form a training set used to train the risk assessment model.

[0092] Model hyperparameter optimization employs the Optuna framework or Bayesian optimization methods for automated search. Searchable hyperparameters include spatial neighborhood window size, number of CNN convolutional layers, number of convolutional kernels, number of LSTM hidden units, number of LSTM layers, number of neurons in the fusion layer, Dropout ratio, learning rate, batch size, loss function type, and optimizer type. Training epochs for each hyperparameter combination are set to 50 to 150 epochs. An early stopping mechanism is triggered when the validation set loss no longer decreases for 5 to 30 consecutive epochs, halting training for the current hyperparameter combination. The gradient clipping threshold is set to 5 to prevent gradient explosion. Optimization objectives can be selected from validation set AUC, F1 score, recall, precision, or a combined score, with recall and F1 score being the primary reference metrics for high-risk samples. The Bayesian optimization method uses a Gaussian process as a surrogate model, selecting the next hyperparameter combination through a sampling function, balancing exploration and utilization.

[0093] Furthermore, the visualization described in step 6 includes:

[0094] Risk probability grid map, risk level zoning map, list of key risk grids, risk factor contribution results, time series risk change curves and model uncertainty results;

[0095] The risk assessment results will be integrated into the geographic information system platform and the 3D visualization platform for comprehensive display.

[0096] Furthermore, the method also includes performing uncertainty quantification analysis on key monitoring points and evaluating the stability of the points by calculating the standard deviation of time-series fluctuations;

[0097] When the standard deviation of the time series fluctuation is greater than a preset threshold, the region is marked as a high uncertainty region and a prompt is made to prioritize manual on-site verification.

[0098] In this embodiment, step 6 involves using the trained CNN-LSTM model for global grid prediction of the study area to obtain the ground collapse risk probability of each grid cell within the target time period. The prediction input includes a static factor feature matrix and a dynamic factor temporal feature matrix, with the matrix dimensions consistent with the input dimensions used during model training. The model's forward inference process employs a batch prediction method, with a batch size set between 64 and 256 to improve inference efficiency. The prediction results are output in the form of CellID and risk probability, where the risk probability is a continuous value between 0 and 1.

[0099] The risk level classification employs a threshold determination method based on the optimal F1 score of the validation set. By traversing different risk probability thresholds and calculating the corresponding F1 scores, the threshold maximizing the F1 score is selected as the optimal threshold. The classification divides the study area into four risk level zones: areas with a risk probability less than 0.25 are classified as low-risk zones; areas with a risk probability between 0.25 and 0.50 are classified as medium-risk zones; areas with a risk probability between 0.50 and 0.75 are classified as relatively high-risk zones; and areas with a risk probability greater than 0.75 are classified as high-risk zones. The thresholds can also be adaptively adjusted based on the actual disaster distribution in the study area. Adjustment methods include quantile thresholding, natural breakpoint methods, and expert experience methods.

[0100] The output includes risk probability raster plots, risk level zoning maps, a list of key risk grids, risk factor contribution results, time-series risk change curves, and model uncertainty results. The risk probability raster plot is output in GeoTIFF format, with grid values ​​ranging from 0 to 1. The risk level zoning map is also output in GeoTIFF format, with grid values ​​ranging from 1 to 4 representing risk level codes. The list of key risk grids is output in CSV format, including CellID, longitude, latitude, risk probability, and risk level fields. Risk factor contribution results are obtained through feature importance analysis methods, including SHAP value analysis or gradient-weighted class activation mapping, to identify the factors with the greatest impact on the risk assessment results. The time-series risk change curve reflects the trend of risk probability changes at the same location at different time points, with the curve data stored in time series format. Model uncertainty results are obtained through Monte Carlo Dropout or ensemble learning methods, outputting the risk probability range with a 95% confidence interval.

[0101] Spatial consistency verification of the prediction results involves comparing and analyzing the predicted results with historical subsidence points, measured deformation areas, and field survey results. Verification metrics include precision, recall, F1 score, AUC, Top-K recall, spatial overlap, and uncertainty interval. Precision is calculated by dividing the number of correctly predicted samples by the total number of samples. Recall is calculated as the proportion of actual high-risk samples identified as high-risk by the model. The F1 score is the harmonic mean of precision and recall. AUC is the area under the ROC curve, reflecting the model's ability to distinguish between positive and negative samples. Top-K recall is calculated as the number of high-risk samples among the K samples with the highest risk probability. Spatial overlap is calculated by dividing the area of ​​spatial overlap between the predicted high-risk area and the actual disaster area by the union area of ​​the two. The spatial consistency error between the evaluation results and the actual situation is controlled within 10%.

[0102] The visualization of risk assessment results is achieved through a 2D GIS platform and a 3D visualization platform. The 2D GIS platform uses software such as ArcGIS, QGIS, or MapGIS to overlay and display risk zoning layers, historical disaster point layers, InSAR deformation layers, rainfall distribution layers, groundwater level change layers, borehole distribution layers, and stratigraphic structure layers. The 3D visualization platform uses software such as Skyline or Cesium to achieve a unified 3D display of surface topography and underground geological structures. Users can use the platform's query function to perform combined queries based on conditions such as regional scope, time period, risk level, and influencing factors. Query results are displayed in the form of tables, charts, or maps. The automatic sorting function for key risk areas sorts grid cells according to risk probability and uncertainty results, prioritizing the display of areas requiring special attention.

[0103] The uncertainty quantification analysis function assesses the risk confidence level of key monitoring points. The assessment method uses the standard deviation of time-series fluctuations, calculated as the standard deviation of the predicted risk probability within a time window. Points with a standard deviation of time-series fluctuations greater than 0.25 are marked as high-uncertainty areas. The system indicates that the geological activity in these areas is relatively intense or the model prediction confidence is low, and recommends prioritizing manual on-site verification. The uncertainty analysis results are published in thematic maps, with red indicating high-uncertainty areas and green indicating low-uncertainty areas.

[0104] Based on the specific implementation of the above steps, the second embodiment of this invention uses the main urban area of ​​a city as the research area to carry out a ground subsidence risk assessment application. The research area covers approximately 500 square kilometers, and the grid division adopts a standard grid of 30 meters by 30 meters, resulting in approximately 550,000 grid units. For data acquisition, the following data were collected: topographic and geomorphological data (DEM data, 30-meter resolution); stratigraphic and lithological data (1:50,000 geological map); geological structural data (fracture distribution map); borehole data (126 borehole points); hydrogeological data (45 groundwater depth monitoring wells); historical subsidence disaster data (62 disaster points); InSAR deformation data (Sentinel-1 data from 2019 to 2023, 96 periods); rainfall data (12 meteorological stations); daily rainfall data; and groundwater level monitoring data (45 monitoring wells, daily monitoring data).

[0105] In terms of data processing, all types of data were unified to the WGS-84 coordinate system and the UTM 51N projection system, resampled to a 30m x 30m grid, and the time series was unified to a 12-day scale. Borehole data underwent spatial interpolation using the ordinary kriging method to generate continuous parameter fields such as water level depth, weak layer thickness, and sand layer thickness, achieving an interpolation accuracy of 85.40%. Outliers and missing values ​​were processed using quality control methods, and all types of data met the modeling requirements.

[0106] For sample construction, high-risk samples were identified based on 62 historical collapse disaster sites, with a buffer radius set to 50 meters, resulting in approximately 200 high-risk samples. Approximately 400 low-risk samples were selected from areas far from disaster sites and with low InSAR deformation. The sample set was divided into training, validation, and test sets in a 7:2:1 ratio.

[0107] For model training, a CNN-LSTM dual-branch architecture was adopted. The CNN branch had three convolutional layers with 32, 64, and 128 kernels respectively, and a 5x5 neighborhood window. The LSTM branch had one bidirectional LSTM layer with 64 hidden units. The fusion layer used an attention fusion mechanism. The initial model used a weighted cross-entropy loss function. After training, the R² was negative, and the prediction accuracy was below 50%, indicating that the model was underfitting and had an imbalanced sample problem.

[0108] For sample augmentation, a GAN was used to generate synthetic high-risk samples. The generator's input noise dimension was 64, and the training rounds were 1000. After screening based on feature distribution similarity, consistency of temporal change trends, and rationality of geological mechanisms, approximately 500 synthetic samples were added, increasing the total number of training samples to approximately 650.

[0109] For hyperparameter optimization, the Optuna Bayesian optimization method was used to search for the optimal combination of hyperparameters. The search space included learning rate, batch size, number of convolutional kernels, number of LSTM hidden units, and Dropout ratio. The final model configuration was: learning rate 0.0005, batch size 32, number of convolutional kernels 64, number of LSTM hidden units 128, and Dropout ratio 0.2, with 100 optimization epochs. After optimization, the model's R² improved to 0.6745, RMSE decreased to 0.6815 mm, and prediction accuracy reached 86.70%. The AUC reached 0.9193, and the mean squared precision (AP) reached 0.8238.

[0110] Regarding risk zoning, the threshold was optimized with the goal of maximizing the F1-Score, and the optimal threshold was set to 0.45. The risk zoning results showed that low-risk areas accounted for approximately 65% ​​of the area, medium-risk areas accounted for approximately 20%, relatively high-risk areas accounted for approximately 10%, and high-risk areas accounted for approximately 5%. Spatial consistency verification results showed that the model achieved a 66.67% capture rate of historical disaster points, and the spatial consistency error was controlled within 10%.

[0111] In terms of visualization applications, the risk assessment results are integrated into the GIS platform and the 3D visualization platform. The platform supports functions such as querying by region, replaying by time, filtering by risk level, and sorting by factor contribution. Uncertainty analysis of key monitoring points shows that there are 3 high-uncertainty points among the 22 key monitoring points, and the system automatically marks them and suggests on-site verification.

[0112] The above specific embodiments demonstrate that the intelligent assessment method for urban ground subsidence provided by the present invention can effectively integrate multi-source data, unify spatiotemporal benchmarks, construct a multi-dimensional factor database, extract spatial and temporal features, train a CNN-LSTM assessment model, achieve sample balancing and hyperparameter optimization, output risk assessment results and visualize them, significantly improving the accuracy, stability and engineering applicability of urban ground subsidence risk assessment.

[0113] According to a third embodiment of the present invention, this embodiment provides an alternative implementation of a ground subsidence risk assessment method for areas with special geological conditions.

[0114] During the data acquisition phase, additional data were added for karst development areas, including karst distribution data, underground river distribution data, and historical karst collapse data. Karst distribution data was obtained through ground-penetrating radar, borehole drilling, and geophysical methods. The data formats include raster-based karst development zoning maps and vector-based cave distribution point maps. Underground river distribution data was obtained through tracer experiments and groundwater level dynamic analysis. The data includes the location of the main underground river channels, flow rate changes, and connectivity evaluation results.

[0115] During the data processing phase, a stratified processing strategy was adopted to address the spatial heterogeneity of karst data. The degree of karst development was categorized into four levels: strongly developed, moderately developed, weakly developed, and undeveloped. Different interpolation methods and parameter settings were used for different levels of regions. Indicator kriging was used for interpolation in strongly developed regions, while ordinary kriging was used in weakly developed regions. The time span of historical karst collapse data was extended to over 30 years to reflect the long-term evolution of the disaster.

[0116] During the model building phase, to address the sudden nature of karst collapses, emergency response features were added as input. These features include the warning level, the operational status of monitoring and early warning equipment, and the adequacy of emergency response measures. The model structure incorporates a multi-task learning branch, simultaneously outputting both risk probability and warning level as objectives. The loss function employs a weighted summation of multi-task losses, with weight coefficients adaptively determined based on validation set performance.

[0117] In the results output phase, a specific assessment of karst collapse is added, including karst collapse susceptibility zoning, the impact range of underground rivers, and emergency response recommendations. Susceptibility zoning is divided into five levels: non-susceptible, low-susceptible, medium-susceptible, high-susceptible, and extremely high-susceptible. The impact range of underground rivers is determined based on the hydraulic connection analysis between the underground river and the surface. Emergency response recommendations are generated according to risk level and warning level, including patrol frequency, monitoring plans, and emergency response measures.

[0118] According to the fourth embodiment of the present invention, this embodiment provides an alternative implementation of the method for assessing the risk of ground subsidence in the impact zone of underground engineering construction.

[0119] During the data acquisition phase, underground engineering data was added, including subway line distribution, subway station locations, shield tunnel alignment, excavation area of ​​the foundation pit, underground utility tunnel alignment, and pile foundation construction records. Subway-related data was obtained from the urban rail transit management department, including the line centerline coordinates, station retaining structure parameters, and shield tunneling parameters. Foundation pit data was obtained from the construction permit management system, including the foundation pit's plan location, excavation depth, support structure, and groundwater level control measures.

[0120] During the data processing phase, underground engineering data is converted into engineering disturbance intensity factors. The calculation method for these factors involves determining dynamic disturbance weights based on the project type, construction stage, and construction progress. The disturbance weight for the shield tunneling stage is set to 0.8, for the foundation pit excavation stage to 0.9, and for the pile foundation construction stage to 0.7. Spatially, the disturbance intensity factors are buffered with the project's influence range as the radius. The influence radius for the shield tunnel section is set to 30 meters, and the influence radius for the foundation pit is set to once the slope height.

[0121] During the sample construction phase, grid cells within the impact area of ​​underground engineering construction were designated as key areas of focus. Sample labels during construction were determined based on on-site monitoring data; samples with deformation rates exceeding warning thresholds were marked as high-risk, while those with deformation rates within control limits were marked as low-risk. Sample labels after construction were determined based on long-term monitoring results, with a follow-up monitoring period of no less than 12 months.

[0122] During the model inference phase, time-varying parameters are added to account for the dynamic changes in engineering construction. These time-varying parameters include construction progress, excavation depth, groundwater level, and monitored deformation rate. The model employs an incremental learning approach, periodically updating model parameters based on newly collected data. Each incremental update involves training for only 10 to 20 epochs to adapt to the dynamic changes in risks during engineering construction.

[0123] In the results output phase, a special risk assessment result for the construction period is added, including construction period risk early warning, settlement development trend prediction, and construction optimization suggestions. The construction period risk early warning uses a four-level system (red, orange, yellow, blue) corresponding to emergency warning, warning, alert, and normal status, respectively. The settlement development trend prediction is based on LSTM time-series prediction capabilities, forecasting settlement changes over the next 1 to 3 months. Construction optimization suggestions are generated based on the risk factor contribution analysis results, including suggestions for adjusting construction parameters and optimizing monitoring schemes.

[0124] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces, or indirect coupling or communication connection between apparatuses or units, and may be electrical, mechanical, or other forms.

[0125] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated units described above can be implemented in hardware or as software functional units. The above are merely embodiments of this application and do not limit the patent scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

[0126] The specific embodiments of the invention have been described in detail above, but they are only examples, and this application is not limited to the specific embodiments described above. For those skilled in the art, any equivalent modifications or substitutions to the invention are also within the scope of this application. Therefore, all equivalent changes, modifications, and improvements made without departing from the spirit and principles of this application should be covered within the scope of this application.

Claims

1. A method for assessing ground subsidence risk through multi-source information fusion, characterized in that, Includes the following steps: Step 1: Collect multi-source data required for urban ground subsidence risk assessment. The multi-source data includes static geological environment data and dynamic environmental disturbance data. Step 2: Perform unified spatiotemporal benchmark transformation and standardization on the multi-source data to unify it to a unified spatial grid; Step 3: Construct a multi-dimensional risk assessment factor database for ground subsidence based on the unified spatial grid; Step 4: Extract spatial features from static factors and extract temporal features from dynamic factors from the risk assessment factor database to construct a risk assessment sample set; Step 5: Construct a CNN-LSTM ground subsidence risk assessment model based on the risk assessment sample set. The model includes a spatial feature extraction branch, a temporal feature extraction branch, a feature fusion layer, and a risk output layer. Step 6: The CNN-LSTM ground collapse risk assessment model is optimized by sample balancing and model hyperparameter optimization, and the risk assessment results are output and visualized.

2. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, The static geological environment data mentioned in step 1 includes: Topographic and geomorphological data, stratigraphic and lithological data, geological structure data, borehole data, hydrogeological data, surface cover data, underground space development data, engineering survey data, road and pipeline network data, and historical ground subsidence disaster data; The dynamic environmental disturbance data includes InSAR surface deformation monitoring data, rainfall data, groundwater level monitoring data, engineering construction disturbance data, surface subsidence rate data, and cumulative subsidence data.

3. The method for assessing ground subsidence risk through multi-source information fusion according to claim 1, characterized in that, Step 2 involves performing unified coordinate transformation, projection processing, cropping, resampling, outlier removal, missing value repair, and field standardization on the multi-source data; unifying all types of data to the same spatial reference system and a unified time scale.

4. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, The database mentioned in step 3 includes: Static factor layer, dynamic factor layer, sample label layer, and model output layer; The static factor layer uses a standard grid as the basic unit to store topographic and geomorphological factors, stratigraphic and lithological factors, hydrogeological factors, geological structural factors, borehole interpolation factors, land cover factors, historical disaster proximity factors, and engineering disturbance background factors. The dynamic factor layer uses grid number and time number as a joint index to store cumulative surface subsidence, subsidence rate over time, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change, and rainfall lag characteristics.

5. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, In step 4, when extracting spatial features from static factors, the values ​​of similar factors within the surrounding neighborhood of the target grid are extracted as spatial neighborhood features. When extracting time-series features of dynamic factors, a time window is constructed with each grid cell as the object to extract cumulative settlement, average settlement rate, maximum settlement rate, cumulative rainfall, maximum daily rainfall, previous rainfall index, groundwater level change amplitude, and rainfall lag term.

6. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, The spatial feature extraction branch in step 5 uses a convolutional neural network to extract the spatial structure features of the target grid and its neighborhood. The time feature extraction branch uses a long short-term memory network to extract the time evolution features of dynamic factors; The feature fusion layer fuses spatial and temporal features; the risk output layer outputs the probability or level of ground subsidence risk.

7. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, The sample balancing process described in step 6 employs at least one of the following methods: undersampling, oversampling, class weight adjustment, and generative adversarial network enhancement. The hyperparameter optimization employs an automated search method to optimize the spatial neighborhood window size, the number of convolutional layers, the number of convolutional kernels, the number of LSTM hidden units, the learning rate, and the batch size.

8. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, In step 6, areas with a risk probability less than the first preset value are classified as low-risk areas, areas with a risk probability greater than the first preset value but less than the second preset value are classified as medium-risk areas, areas with a risk probability greater than the second preset value but less than the third preset value are classified as relatively high-risk areas, and areas with a risk probability greater than the third preset value are classified as high-risk areas.

9. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, The visualization described in step 6 includes: Risk probability grid map, risk level zoning map, list of key risk grids, risk factor contribution results, time series risk change curves and model uncertainty results; The risk assessment results will be integrated into the geographic information system platform and the 3D visualization platform for comprehensive display.

10. The method for assessing ground subsidence risk by multi-source information fusion according to claim 1, characterized in that, It also includes conducting uncertainty quantification analysis on key monitoring points and evaluating the stability of the points by calculating the standard deviation of time-series fluctuations; When the standard deviation of the time series fluctuation is greater than a preset threshold, the region is marked as a high uncertainty region and a prompt is made to prioritize manual on-site verification.

Images