Loess area urban geological disaster risk evaluation method based on underground water level

Abstract

The invention relates to the technical field of geological disaster risk assessment, and discloses a loess area urban geological disaster risk assessment method based on an underground water level. The method comprises the following steps: acquiring underground water level data by arranging a monitoring network, and constructing a space-time sequence after processing; further extracting an interaction mode between the underground water level and the geological parameters to form a coupling feature vector; the method comprises the following steps: dividing a region into sub-blocks of different levels according to data resolution differences, customizing grid specifications, establishing a multi-resolution grid framework, realizing space distribution and time expansion of feature vectors, and generating a multi-resolution feature field; extracting risk indicators and combining the risk indicators into a risk spectrum; and finally, solving risk values of the key points through a risk calculation engine, and matching risk levels in combination with a knowledge base to form a regional evaluation conclusion. According to the invention, through coupling feature extraction and a multi-resolution grid adaptive technology, the accuracy and reliability of risk evaluation are improved.

Description

Technical Field

[0001] This invention relates to the field of geological disaster risk assessment technology, specifically a method for assessing urban geological disaster risks in loess areas based on groundwater levels. Background Technology

[0002] The occurrence of urban geological disasters in loess regions is closely related to groundwater level fluctuations. Existing risk assessment methods typically employ a uniform grid scale for spatial analysis. This homogeneous grid division method is ill-suited to the reality of uneven distribution of monitoring points and significant differences in data resolution within urban areas. In data-dense areas, fixed large grids lose a significant amount of detailed information, failing to accurately capture local high-risk points; while in data-sparse areas, overly dense grids lead to insufficient data support, introducing substantial uncertainty or computational redundancy. This rigidity in spatial processing is a major drawback of existing technologies.

[0003] Traditional methods for analyzing the impact of groundwater levels on geological hazards often treat groundwater level parameters as relatively independent or simply linearly superimposed with geological parameters such as soil properties and stratigraphic structure. This approach fails to reveal the complex nonlinear interaction mechanism between dynamic fluctuations in groundwater levels and the response of geological bodies. Due to the lack of a quantitative characterization of this intrinsic coupling relationship, the model struggles to accurately reflect the instability process of geological bodies under the influence of groundwater, resulting in insufficient explanation of the risk formation mechanism and affecting the accuracy of the evaluation results.

[0004] A technique is needed that can adapt to spatial heterogeneity of data and perform multi-source parameter coupling analysis to solve the above problems and achieve accurate quantitative evaluation of geological disaster risks in Loess urban areas. Summary of the Invention

[0005] The purpose of this invention is to provide a method for assessing urban geological disaster risks in loess areas based on groundwater levels, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, this invention provides a method for assessing urban geological hazard risk in loess areas based on groundwater level, the method comprising: A monitoring equipment network was deployed to continuously collect raw groundwater level data in urban areas of the Loess Plateau region; noise filtering and missing value imputation were performed on the raw groundwater level data to generate a clean monitoring dataset. By processing the clean monitoring dataset, an initial spatiotemporal sequence with a unified format is constructed through time alignment and spatial registration. Feature extraction techniques are then used to analyze the interaction pattern between groundwater level and geological parameters from the initial spatiotemporal sequence, and a coupled feature vector is output. The investigation examines the resolution distribution of the monitoring data and divides the region into sub-blocks of different levels based on the resolution. A customized grid specification, including cell size and grid density, is then established for each sub-block. Initiate the mesh generation process, create a multi-resolution mesh framework based on customized mesh specifications, and perform spatial allocation and temporal expansion of coupled feature vectors within the multi-resolution mesh framework to generate a multi-resolution feature field. Examine the time series of multi-resolution feature fields, extract risk indicators using pattern recognition methods, and combine them into a risk spectrum; The geological hazard risk calculation engine is invoked, and the risk spectrum is used as input to solve for the risk value of key geological hazard points, obtain the risk status dataset, scan the risk status dataset, detect risk exceeding the limit events of each key geological hazard point, and record the number of events. Access the geological disaster risk knowledge base, match risk curves based on the frequency of events, assign risk levels, integrate the risk levels of each key point, and form a regional risk assessment conclusion.

[0007] Preferably, the step of performing noise filtering and missing value imputation on the raw groundwater level data to generate a clean monitoring dataset includes: The system receives raw groundwater level data streams, uses a sliding window algorithm to detect abnormal fluctuations, removes these abnormal fluctuations, and uses interpolation to fill in the data gaps. The processed data is standardized to adjust the data scale and distribution, making the data from all monitoring points comparable. Verify the integrity and consistency of the data, mark invalid data segments; remove invalid data segments, and finally output the clean monitoring dataset.

[0008] Preferably, the process of processing the clean monitoring dataset involves constructing an initial spatiotemporal sequence with a uniform format through time alignment and spatial registration, including: Read the cleanroom monitoring dataset, unify the time sampling frequency, and resample data from different frequencies to a common time axis; Align the spatial coordinate system to map the monitoring point locations to a standard geographic grid; extract the groundwater level time series for each grid point; Aggregate the time series of all grid points and organize them into an initial spatiotemporal sequence matrix according to time order and spatial location.

[0009] Preferably, the step of using feature extraction technology to analyze the interaction pattern between groundwater level and geological parameters from the initial spatiotemporal sequence and outputting a coupled feature vector includes: Load the initial spatiotemporal sequence, combine it with the geological parameter database, including soil type and slope information; calculate the joint statistics of groundwater level and geological parameters; Principal component analysis is used to reduce the feature dimensionality while preserving the main direction of variation. Time sliding window analysis is then performed on the dimensionality-reduced features to enhance dynamic interaction signals. The enhanced feature set is then compiled to generate coupled feature vectors.

[0010] Preferably, the resolution distribution of the survey and monitoring data is such that the region is divided into sub-blocks of different levels according to the resolution; and a customized grid specification is provided for each sub-block, including cell size and grid density, including: Analyze the spatial resolution of the monitoring data and draw resolution contour maps; divide the region into high-resolution sub-blocks, medium-resolution sub-blocks, and low-resolution sub-blocks based on the contour maps; Set small cell size and high grid density for high-resolution sub-blocks; set medium cell size and medium grid density for medium-resolution sub-blocks; set large cell size and low grid density for low-resolution sub-blocks. Generate a grid specification configuration file for each sub-block.

[0011] Preferably, the initiation of the mesh generation process, which creates a multi-resolution mesh framework based on customized mesh specifications, includes: Parse the mesh specification configuration file and initialize the mesh generator; create mesh patches independently for each sub-block and record the patch vertices and connection relationships; Assemble the facets into a 3D mesh and assign volume element identifiers; check the continuity and boundary consistency of the mesh. Optimize the mesh structure, smooth the mesh transition areas, and output a multi-resolution mesh framework.

[0012] Preferably, the step of performing spatial allocation and temporal expansion of coupled feature vectors within a multi-resolution mesh framework to generate a multi-resolution feature field includes: Project the coupled feature vectors onto the vertices of the multi-resolution mesh frame; fill the feature values ​​within the mesh cells using a spatial interpolation algorithm; Extend the feature values ​​along the time axis to generate a feature snapshot at each time step; serialize all feature snapshots to construct a time series of a multi-resolution feature field.

[0013] Preferably, the step of examining the time series of multi-resolution feature fields, extracting risk indicators using pattern recognition methods, and combining them into a risk spectrum includes: Decompose the time series of the multi-resolution feature field into trend components, periodic components, and residual components; calculate the statistical properties of each component; Identify patterns related to geological hazards, such as abrupt change points and periodic oscillations, quantify the intensity of these patterns, form a set of risk indicators, arrange the risk indicator set in chronological order, and synthesize a risk spectrum.

[0014] Preferably, the step of calling the geological hazard risk calculation engine, taking the risk spectrum as input, and solving for the risk values ​​of key geological hazard points to obtain a risk status dataset includes: Configure the parameters of the geological hazard risk calculation engine, including material properties and boundary conditions; input the risk spectrum as an external stimulus; Perform numerical simulations to calculate the stress, strain, or displacement response at each key point of the geological hazard; extract risk indicators from the response data. Compile risk indicators for all key points to generate a risk status dataset.

[0015] Preferably, the scanned risk status dataset detects risk exceeding events at each key geological hazard point and records the number of events occurring; it accesses the geological hazard risk knowledge base, matches risk curves based on the number of events occurring, and assigns risk levels, including: Set a risk threshold, iterate through the risk status dataset, and mark risk values ​​that exceed the threshold; count the number of times each key point is marked as the number of events. Query the geological disaster risk knowledge base and retrieve the preset risk curve; input the number of times the event occurred, read the corresponding risk level from the risk curve, assign a risk level to each key point, and store the evaluation results.

[0016] Compared with the prior art, the beneficial effects of the present invention are: By investigating the resolution distribution of monitoring data, the region is hierarchically divided, and differentiated grid specifications are customized for different sub-blocks, generating a multi-resolution grid framework. This method abandons the constraints of traditional fixed grids, achieving adaptive matching between grid scale and data density. Small-size, high-density grids are used in high-resolution data areas to fully preserve local details; large-size, low-density grids are used in low-resolution data areas to avoid excessive interpolation without data support, effectively balancing computational efficiency and model accuracy. This data-driven gridding approach makes the spatial analysis process more closely resemble actual data conditions, improves the rationality of spatial representation, and reduces errors caused by inappropriate grid scale.

[0017] Feature extraction techniques are employed to analyze the interaction patterns of groundwater levels and geological parameters in a spatiotemporal sequence, outputting coupled feature vectors. This technique transcends simple correlations between parameters, capturing dynamic correlations such as the relationship between groundwater level rise / fall rates and soil strength parameter decay, and pore water pressure dissipation and formation stability. This deep feature coupling integrates previously isolated environmental factors into comprehensive indicators reflecting the system's intrinsic behavioral mechanisms, enabling risk assessment models to more fundamentally characterize the gestation process of geological hazards. The resulting risk spectrum focuses more on revealing the hazard-causing mechanisms, enhancing the theoretical basis and physical significance of the assessment results, and improving the reliability of risk state identification. Attached Figure Description

[0018] Figure 1 This is a schematic diagram illustrating the working principle of the urban geological hazard risk assessment method based on groundwater level in the loess region described in this invention. Figure 2 A flowchart for noise filtering and missing value imputation; Figure 3 The flowchart for extracting coupled feature vectors; Figure 4 Scatter plot showing the correlation between the principal component distribution and risk contribution of geological hazard coupling characteristics in the Loess region; Figure 5 A bar chart showing the risk level distribution of key geological disaster points in the Loess Plateau region. Detailed Implementation

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

[0020] Please see Figure 1 This invention provides a method for assessing urban geological hazard risk in loess areas based on groundwater level. The method includes: deploying a monitoring equipment network to continuously collect raw groundwater level data in urban areas of loess regions. The monitoring equipment includes pressure level gauges, conductivity sensors, etc., and records the groundwater level depth at minute or hourly intervals. The raw groundwater level data is then subjected to noise filtering and missing value imputation to generate a clean monitoring dataset. A sliding window algorithm is used to identify abnormal fluctuations caused by equipment failure or environmental interference, and linear interpolation or spline interpolation is used to fill data gaps, ensuring the continuity of the time series. The clean monitoring dataset is processed by time alignment and spatial registration to construct an initial spatiotemporal sequence with a unified format, unifying data from different monitoring points to the same timestamp and geographic coordinate system. Feature extraction technology is used to analyze the interaction patterns between groundwater level and geological parameters from the initial spatiotemporal sequence, outputting a coupled feature vector. Geological parameters include loess porosity, permeability coefficient, and cohesion.

[0021] The investigation examines the resolution distribution of monitoring data, dividing the region into sub-blocks of varying resolution levels, which depend on the monitoring point density and data acquisition frequency. A customized grid specification, including cell size and grid density, is established for each sub-block; a fine grid is used for high-resolution areas, while a sparse grid is used for low-resolution areas. A grid generation process is initiated to create a multi-resolution grid framework based on the customized grid specifications. The grid type can be unstructured or structured. Within the multi-resolution grid framework, spatial allocation and temporal expansion of coupled feature vectors are performed to generate a multi-resolution feature field. Spatial interpolation is performed using Kriging interpolation or inverse distance weighting. The time series of the multi-resolution feature field is examined, and risk indicators are extracted using pattern recognition methods and combined to form a risk spectrum. Risk indicators include the rate of change of groundwater level and the amplitude of periodic fluctuations. A geological hazard risk calculation engine is invoked, using the risk spectrum as input, to calculate the risk values ​​of key geological hazard points, obtaining a risk status dataset. Key points include slope toes and the area surrounding building foundation pits. The system scans the risk status dataset to detect risk exceedance events at each key geological hazard point, recording the frequency of these events. Exceedance events are defined as risk values ​​exceeding preset safety thresholds. It accesses the geological hazard risk knowledge base, matches risk curves based on event frequency, and assigns risk levels. These risk curves are built based on historical hazard cases and numerical simulation results. Finally, it integrates the risk levels of each key point to form a regional risk assessment conclusion, outputting a risk level distribution map or risk assessment report.

[0022] Example 1: See Figure 2The system receives raw groundwater level data streams from a network of monitoring equipment deployed in urban areas of loess regions. The data format includes a timestamp field, a unique identifier for each monitoring point, and a numerical field for groundwater level elevation or depth. In implementation, a sliding window algorithm is used to detect abnormal fluctuations. This algorithm moves across the data stream at fixed time intervals, configurable to 24 or 48 hours. The algorithm calculates the arithmetic mean and standard deviation of the data within each window, identifying data points deviating from the arithmetic mean by more than three times the standard deviation as abnormal fluctuations. Abnormal fluctuations are removed, and interpolation methods are used to fill data gaps. Cubic spline interpolation is preferred to ensure curve smoothness, but linear interpolation can also be used in areas with large data gaps or boundary regions. The processed data undergoes standardization using the Z-score method, which involves subtracting the arithmetic mean of the entire monitoring point's time series from each data point and then dividing by the standard deviation of that time series. This adjusts the data scale and distribution, ensuring that all monitoring points have zero mean and unit variance, thus achieving direct comparability. Verify the integrity and consistency of the data, checking for consecutive missing data points in the time series, whether the values ​​exceed physically reasonable ranges, and whether there are contradictions in overlapping periods between different monitoring points. Mark invalid data segments, which are data intervals with consecutive missing data points exceeding a preset number or whose values ​​clearly violate hydrogeological patterns. Remove invalid data segments, ultimately outputting a clean monitoring dataset. The clean monitoring dataset is stored in a structured table or time series database format, ensuring that each valid data point has an accurate timestamp, location information, and corrected water level value.

[0023] The cleanroom monitoring dataset is read and serves as the foundational data source for subsequent processing steps. A unified time sampling frequency is established, resampling data from different frequencies onto a common time axis. In practice, this common time axis is set to one sampling point per hour. Data with sampling frequencies higher than hourly undergo downsampling, for example, aggregating minute-level data into hourly data using arithmetic mean or median filtering. Data with sampling frequencies lower than hourly undergo upsampling, for example, interpolating daily data into hourly data using linear interpolation or spline interpolation. The spatial coordinate system is aligned, mapping monitoring point locations to a standard geographic grid. The spatial coordinate system uniformly adopts the national geodetic coordinate system, such as CGCS2000. The standard geographic grid uses a regular raster, with the grid size determined based on the area and data density. The groundwater level time series for each grid point is extracted. For grids with monitoring points, the resampled time series for that monitoring point is directly used. For grids without monitoring points, inverse distance weighted interpolation or kriging interpolation is used to generate a virtual time series based on data from surrounding monitoring points. The time series data of all grid points are aggregated and organized into an initial spatiotemporal sequence matrix according to chronological order and spatial location. This initial spatiotemporal sequence matrix is ​​a three-dimensional array. The first dimension represents the time point, and the second and third dimensions together represent the planar coordinates of the grid points. Array elements store the groundwater level values ​​at the corresponding spatiotemporal locations. In some embodiments, the initial spatiotemporal sequence matrix can be stored in scientific data formats such as NetCDF or HDF5 to support efficient multidimensional data reading, writing, and manipulation. It can be understood that the construction of the initial spatiotemporal sequence matrix provides a structured data foundation for subsequent feature extraction and spatiotemporal analysis.

[0024] In the data preprocessing stage, the parameter settings of the sliding window algorithm have a significant impact on the sensitivity of anomaly detection. A window length that is too short may lead to misjudgment of normal fluctuations, while a window that is too long may mask true anomalies. In practice, the window length needs to be adjusted according to the natural cycle of groundwater level changes, such as considering diurnal or seasonal variations. The choice of interpolation method also requires careful consideration. Cubic spline interpolation can better maintain the continuity and smoothness of water level changes, but it may cause overshoot at data gap edges. Linear interpolation, while simple and stable, may introduce broken-line artificial traces. Standardization transformation is a crucial step in ensuring data comparability. The Z-score method requires data to approximately follow a normal distribution. When the distribution of monitoring point data deviates significantly from normality, it may be necessary to consider using quartile ranges for standardization or performing data transformation. When verifying data integrity and consistency, in addition to automated rule checks, manual verification should be conducted in conjunction with hydrogeological expertise. For example, it should be determined whether sudden changes in water level are consistent with known precipitation events or artificial pumping activities to avoid mistakenly deleting genuine but reasonable anomalies.

[0025] Resampling can introduce errors when constructing the initial spatiotemporal series. Aggregation methods used in downsampling can affect data representativeness; arithmetic mean reflects the overall level but may smooth out peaks, while median filtering is insensitive to outliers but may lose some detailed information. Interpolation in upsampling is essentially a prediction based on assumptions, and its accuracy depends on the continuity and stability of the groundwater system. When spatially mapping to a standard geographic grid, the choice of grid size requires a trade-off between resolution and computational cost; too small a grid will result in many grids having no data for interpolation, while too large a grid will lose spatial details. The quality of the generated virtual time series is highly dependent on the choice of spatial interpolation algorithm and parameter settings. Inverse distance weighting is simple to implement but may produce a "bull's-eye" effect, while Kriging interpolation provides the best unbiased estimate but requires calculating a variogram model. The organization of the initial spatiotemporal series matrix directly affects the efficiency of subsequent algorithms; time-priority storage facilitates time series analysis, while space-priority storage facilitates spatial analysis. A suitable data layout needs to be selected based on the main computational task. Understandably, each technology choice needs to be tested and optimized in specific application scenarios to achieve a balance between accuracy and efficiency.

[0026] Optionally, more advanced denoising algorithms, such as wavelet transform thresholding, can be introduced in the data preprocessing stage. Wavelet transform can analyze signals in both the time and frequency domains simultaneously, better separating noise and useful information. For missing value imputation, with sufficient historical data, time series prediction models such as the ARIMA model can be considered for predictive imputation to improve the accuracy of imputation in cases of long-term missing values. When constructing spatiotemporal series, for spatial interpolation, co-kriging interpolation can be performed by combining auxiliary variables such as elevation and soil type to utilize more information and improve interpolation accuracy. The initial spatiotemporal series matrix can also be enhanced with spatiotemporal indexing information, such as recording the geographic coordinates and elevation of each grid point, to facilitate more complex spatiotemporal queries and spatial analyses. These enhancement methods can be selectively implemented based on actual data quality and computational resources.

[0027] Example 2: See Figure 3The initial spatiotemporal sequence is structured data organized in the form of a three-dimensional array, containing time dimension, spatial dimension X, spatial dimension Y, and corresponding groundwater level observations. Combined with a geological parameter database, which is a static dataset obtained through field geological surveys, laboratory geotechnical tests, and remote sensing interpretation, the database contains fields such as soil type code, saturation percentage, internal friction angle, and slope degree. The database is stored in a raster layer or vector layer format supported by the geographic information system, with each raster cell or vector area associated with corresponding geological attribute values. The joint statistics of groundwater level and geological parameters are calculated. In practice, the calculation of joint statistics is carried out within each standard geographic grid cell. For each grid cell, the groundwater level time series and the geological parameter values ​​of its location are extracted. The calculation items include Pearson correlation coefficient, Spearman rank correlation coefficient, covariance matrix, and mutual information entropy based on information theory. Pearson correlation coefficient quantifies the degree of linear correlation between groundwater level and continuous geological parameters. Spearman rank correlation coefficient is used to assess the monotonic relationship with ordered geological parameters. Covariance matrix describes the magnitude and direction of the co-change of multiple geological parameters and groundwater level. Mutual information entropy can capture the mixed linear and nonlinear correlation relationship.

[0028] Principal Component Analysis (PCA) is used to reduce the dimensionality of the features. The input to PCA is an initial feature matrix composed of all the joint statistics calculated above. Rows in the initial feature matrix correspond to different grid cells, and columns correspond to different joint statistic indices. PCA calculates the covariance matrix of the initial feature matrix, and then solves for the eigenvalues ​​and eigenvectors of the covariance matrix. The magnitude of the eigenvalues ​​reflects the amount of variability in the original data contained in each principal component direction. The main variability directions are retained. Specifically, based on the cumulative contribution rate of the eigenvalues, a threshold is set, such as 85% or 90%. Principal components are selected starting from the eigenvalue with the largest contribution rate until the cumulative contribution rate of the selected principal components exceeds the set threshold. These selected principal component directions represent the most significant variation patterns in the original joint statistics. A time-sliding window analysis is then performed on the dimensionality-reduced features. The dimensionality-reduced features, i.e., the scores of each grid cell on the principal components, constitute a new, lower-dimensional time series dataset. Time-sliding window analysis slides along the time axis with a fixed window length, set according to the research objective (e.g., 30 days, 90 days). Within each window, the statistical properties of the principal component scores are calculated, such as the arithmetic mean, variance, slope of the linear trend, and autocorrelation coefficient. These window statistics aim to enhance the capture of dynamic interaction signals between groundwater level and geological parameters. For example, the trend slope can reflect the long-term direction of change in the interaction relationship, and the variance can reflect the stability of the interaction strength. The enhanced feature set is compiled to generate coupled feature vectors. The enhanced feature set contains the principal component scores and their window statistics for each grid cell in each time window. This information is arranged in a fixed order, and each coupled feature vector uniquely corresponds to a comprehensive state description of a grid cell within a specific time window.

[0029] In some embodiments, the spatial heterogeneity and scale effects of geological parameters may need to be considered when calculating joint statistics. For example, for a grid cell containing multiple soil types, the area-weighted average correlation between groundwater level and the main soil type within that grid cell can be calculated, or the maximum value of the correlation statistics with each soil type can be calculated to identify the most sensitive geological factor. Data standardization before principal component analysis is an important step. Since different joint statistics have different dimensions and numerical ranges, they need to be standardized to have zero mean and unit variance to avoid large-scale statistics dominating the results of principal component analysis. The length of the time sliding window needs to be carefully chosen; a window that is too long will smooth out short-term, dramatic interaction signals, while a window that is too short may lead to unstable statistics and excessive noise. It can be understood that the generation of coupled feature vectors is a process of data condensation and feature enhancement, which transforms the original, high-dimensional spatiotemporal sequence data into a lower-dimensional representation that better reveals the underlying mechanisms and is suitable for subsequent risk pattern recognition.

[0030] Optionally, when calculating joint statistics, in addition to the conventional correlation coefficient, more complex measures can be introduced, such as the maximum information coefficient, which can detect a wider range of functional relationships beyond just linear or monotonic relationships. In principal component analysis, besides selecting principal components based on cumulative contribution rates, statistical methods such as scree plots or parallel analysis can be combined to determine the number of principal components to retain, thereby improving the robustness of the results. For time-sliding window analysis, the window function is not limited to rectangular windows; weighted windows such as Hanning windows or Hamming windows can be used to reduce spectral leakage. The window step size can also be set to be smaller than the window length to achieve overlapping sliding, thereby generating a more temporally continuous feature sequence. The compilation of coupled feature vectors can introduce feature cross operations, such as multiplying or dividing the window statistics of different principal components, to explicitly construct interactive features, which may capture more complex nonlinear interaction patterns. These optional methods can be selectively applied according to specific data characteristics and computational resources.

[0031] It is understandable that the quality of the coupled feature vector directly affects the accuracy of subsequent risk assessment. To quantitatively assess the amount of information retained in the feature vector, an information retention index can be calculated. An exemplary formula for calculating the information retention index is as follows:

[0032] in: This represents the degree of information retention; the closer the value is to 1, the more complete the information is retained. This indicates the total number of time windows. This indicates the total number of spatial grid cells. Indicates the first The first time window, the first At each grid cell, a certain comprehensive representation value of the high-dimensional original eigenvector composed of all original joint statistics, such as its magnitude or the score of the first principal component. Indicates the first The first time window, the first At each grid cell, the corresponding comprehensive representation value of the finally generated coupled feature vector. It is a very small positive number used to prevent the denominator from being zero.

[0033] See Figure 4This figure is a visualization of the coupled feature vector analysis process. It presents the correlation between the principal component scores and risk contribution of each grid cell within the study area in scatter plot form: the horizontal axis represents the principal component 1 score, the vertical axis represents the principal component 2 score, and the color maps the risk contribution. The core function of this figure is to intuitively identify high-risk feature clusters, providing a spatial distribution basis for the subsequent extraction of the risk spectrum. It is a key visualization tool for extracting risk-related features from high-dimensional spatiotemporal data.

[0034] Example 3: The specific quantitative indicator of spatial resolution is the number of effective monitoring points per unit area. The calculation method involves using geographic information system (GIS) software to perform kernel density estimation on the monitoring point distribution, generating a continuous density surface raster map. The value of each raster cell represents the monitoring point density within a certain search radius. A resolution contour map is then drawn based on the raster map generated by the kernel density estimation. Contour lines connect areas with the same monitoring point density, forming closed or open curves. The numerical intervals of the contour maps need to be reasonably set according to the actual density distribution range to clearly reflect the resolution change trend. Based on the contour map, the region is divided into high-resolution, medium-resolution, and low-resolution sub-blocks. The division threshold is determined based on the overall density distribution characteristics of the monitoring network. For example, areas with a monitoring point density higher than 10 per square kilometer are classified as high-resolution sub-blocks, areas with a density between 3 and 10 per square kilometer are classified as medium-resolution sub-blocks, and areas with a density lower than 3 per square kilometer are classified as low-resolution sub-blocks. The boundaries of the sub-blocks strictly follow the contour lines.

[0035] For high-resolution sub-blocks, small cell sizes and high grid densities are used. In practice, the grid cell size for high-resolution sub-blocks can be set to 10 meters by 10 meters, meaning each grid cell covers 100 square meters of land area, with a corresponding grid density of 10,000 cells per square kilometer. This fine division aims to fully utilize high-density monitoring data to capture local detail changes. For medium-resolution sub-blocks, medium cell sizes and medium grid densities are used. The grid cell size for medium-resolution sub-blocks can be set to 50 meters by 50 meters, with each cell covering 2,500 square meters, and a grid density of 400 cells per square kilometer. This setting balances computational efficiency and spatial representation accuracy. For low-resolution sub-blocks, large cell sizes and low grid densities are used. The grid cell size for low-resolution sub-blocks can be set to 100 meters by 100 meters, with each cell covering 10,000 square meters, and a grid density of 100 cells per square kilometer. This relatively sparse grid is suitable for data-scarce areas and avoids uncertainties introduced by excessive interpolation. Generate a mesh specification configuration file for each sub-block. The mesh specification configuration file is written in JSON format and contains a unique identifier for the sub-block, the type of the sub-block (high, medium, or low resolution), the coordinate string of the boundary polygon of the sub-block, the specified mesh cell size (length and width), the calculated mesh density, and a description of the mesh type (such as a regular rectangular mesh).

[0036] The process involves parsing the mesh specification configuration file and initializing the mesh generator. The mesh generator is a dedicated software module for creating spatial meshes; it reads parameters from the configuration file to determine the meshing rules for each sub-block. Mesh patches are created independently for each sub-block. For regular rectangular meshes, the mesh patch creation process first determines a minimum bounding rectangle based on the sub-block boundary, then divides the rectangle into regular rows and columns according to the specified cell size. The vertices of the mesh patch are the corner points of each mesh cell, and the connection relationships are naturally defined by the row and column indices of the cells; adjacent cells share vertices. The patches are then assembled into a 3D mesh volume. The 3D mesh volume is constructed by assigning an elevation value to each surface mesh patch. The elevation data comes from the region's digital elevation model, thus elevating the 2D planar mesh into a 3D mesh volume. Each 3D mesh cell is assigned a unique volume cell identifier, typically a long integer. The mesh continuity and boundary consistency are checked. Continuity checks ensure there are no voids or gaps between mesh cells. Boundary consistency checks focus on adjacent areas of different resolution sub-blocks, ensuring seamless connection between the meshes of high-resolution sub-blocks and adjacent low-resolution sub-blocks at the boundaries, with vertices and edges correctly matched without overlap or breaks. The mesh structure is optimized, smoothing transition regions. Transition regions refer to the connecting areas between different resolution sub-blocks. Optimization operations include locally refining or coarsening mesh cells in transition regions, and applying the Laplacian smoothing algorithm to fine-tune the positions of mesh nodes to reduce sharp transitions and make mesh changes smoother. A multi-resolution mesh framework is output and stored in a standard mesh file format.

[0037] In some embodiments, the sub-block division may not be a simple three-level system, but rather more levels may be set according to actual needs, such as adding ultra-high resolution sub-blocks (for key monitoring areas) and extremely low resolution sub-blocks (for areas with extremely sparse data). The shape of the grid cells is not limited to rectangles; in complex boundary areas, triangular or irregular polygonal grid patches can be used to better fit natural boundaries. The construction of the 3D mesh can consider not only surface elevation but also underground geological strata information, such as dividing the mesh into several layers in the vertical direction to construct a true 3D geological model for analyzing the impact of groundwater levels at different depths. When checking the mesh quality, in addition to continuity and boundary consistency, indices such as the Jacobian matrix or aspect ratio of the grid cells can be calculated to evaluate whether the shape of the grid cells is beneficial to subsequent numerical calculations. It can be understood that the core of the multi-resolution mesh framework design is to invest more computational resources in data-intensive areas to obtain accuracy and reduce computational costs in data-sparse areas, thereby achieving an optimized balance between overall computational efficiency and accuracy.

[0038] Optionally, spatial autocorrelation analysis can be introduced as an auxiliary criterion when dividing sub-blocks. For example, Moran's index can be used to assess the spatial clustering of monitoring point data. Areas with high spatial autocorrelation may be divided into higher-resolution sub-blocks even if the monitoring density is slightly lower, thus maintaining the integrity of geological features. The grid generation process can employ adaptive grid densification technology. Based on an initial uniform grid, the local grid density is dynamically adjusted according to the gradient magnitude of the coupled feature vectors, automatically densifying the grid in areas with drastic feature changes. For smoothing grid transition regions, more advanced grid smoothing algorithms, such as angle-weighted or spring-model-based smoothing methods, can be used to obtain better grid quality. The output of the multi-resolution grid framework can include detailed metadata describing the grid generation parameters, coordinate system information, and data source descriptions for each sub-block, facilitating subsequent traceability and reuse. These optional methods can further enhance the scientific rigor and practicality of the grid framework.

[0039] To quantitatively evaluate the rationality of sub-block partitioning schemes with different resolutions, an evaluation function for the partitioning scheme can be defined. An exemplary evaluation function for a partitioning scheme is as follows:

[0040] in: The smaller the absolute value, the better the partitioning scheme reflects the natural distribution of data density. This represents the total area of ​​the entire study region. This indicates the total number of sub-blocks that have been divided. Indicates the first The area of ​​each sub-block. Indicates the first Average monitoring point density within each sub-block. This represents the average density of monitoring points across the entire study area. Indicates the first The variance of the density of monitoring points within each sub-block, This represents the average value of the global density variance.

[0041] Example 4: Projecting the coupled feature vectors onto the vertices of the multi-resolution mesh framework. The coupled feature vectors are the feature set generated in the previous steps, containing principal component scores and window statistics. The vertices of the multi-resolution mesh framework are the basic spatial nodes constituting the mesh. The projection operation is based on the spatial nearest neighbor matching principle. For each vertex in the multi-resolution mesh framework, the spatial Euclidean distance to all monitoring points (i.e., the source locations of the coupled feature vectors) is calculated, and the coupled feature vector corresponding to the nearest monitoring point is assigned to that mesh vertex. Spatial interpolation algorithms are used to fill the feature values ​​within the mesh cells. In specific implementations, for each mesh cell (such as a tetrahedron or hexahedron) of the multi-resolution mesh framework, the feature values ​​of its internal points need to be obtained by interpolation using the known feature values ​​of its vertices. The commonly used spatial interpolation algorithm is linear interpolation. For hexahedral cells, it is assumed that the feature values ​​within the cell change linearly in the three coordinate axes. The feature values ​​of any point within the cell are calculated using the trilinear interpolation formula based on the feature values ​​of the eight vertices. The eigenvalues ​​are extended along the time axis to generate a feature snapshot at each time step. Each snapshot represents the complete state of the entire study area at a specific moment, containing the interpolation results of all feature components for each grid cell at that moment. All feature snapshots are serialized to construct a time series of the multi-resolution feature field. Serialization refers to storing the feature snapshots of each time step sequentially in a large multidimensional array or a dedicated time series database, forming a multi-resolution feature field data cube that is continuous in time and space.

[0042] The time series of a multi-resolution feature field is decomposed into trend, periodic, and residual components. The time series of the multi-resolution feature field is a multivariate time series for each grid cell. The decomposition process uses a seasonal decomposition algorithm, such as STL decomposition. STL decomposition iteratively decomposes the time series into trend, periodic, and residual components using a locally weighted regression scatter smoothing method. The trend component captures the long-term, slow direction of change of the features; the periodic component extracts fluctuation patterns with fixed periods; and the residual component contains the random fluctuations and noise remaining after removing the trend and periodicity. The statistical properties of each component are calculated. For the trend component, its linear regression slope is calculated to quantify the rate of change, and its fluctuation range is calculated to quantify the amplitude of change. For the periodic component, its amplitude is calculated to quantify the intensity of fluctuation, and its dominant frequency is calculated to identify the dominant period. For the residual component, its standard deviation is calculated to quantify the magnitude of random fluctuations, and its skewness and kurtosis are calculated to understand the distribution shape. Patterns related to geological hazards, such as abrupt changes and periodic oscillations, are identified. Abrupt change detection employs a cumulative sum algorithm or a Bayesian change point detection algorithm to identify moments in the feature sequence where statistical characteristics change significantly. Periodic oscillation analysis uses a Fast Fourier Transform to convert periodic components to the frequency domain, identifying frequency peaks with energy significantly higher than background noise; the periods corresponding to these peaks are important oscillation modes. The intensity of these modes is quantified to form a set of risk indicators. The intensity of an abrupt change is represented by the difference in the mean or variance of the sequence before and after quantification. The intensity of a periodic oscillation is represented by the ratio of its amplitude to the amplitude of the background noise. These quantified indicators, along with their location, time, and corresponding feature component information, are encapsulated into structured risk indicators. The risk indicator set is arranged chronologically to synthesize a risk spectrum. The risk spectrum is a multidimensional time series, with each time point corresponding to a vector containing the type, spatial location, and intensity information of all activated risk indicators at that moment. The risk spectrum is stored in tabular or timestamped key-value pair format.

[0043] In some embodiments, spatial interpolation at the boundaries of sub-blocks with different resolutions within a multi-resolution grid framework requires special handling to ensure a smooth transition of the feature field across resolution boundaries. Constrained interpolation methods can be employed, introducing continuity conditions at the boundaries to prevent jumps in feature values. When the time series of the multi-resolution feature field contains a large number of missing values, traditional methods such as STL decomposition may fail. Robust decomposition algorithms capable of handling missing values ​​are needed, or the missing values ​​can be reasonably imputed before decomposition. The quantification of risk indicators may not rely solely on a single feature component but requires the integration of pattern information from multiple feature components. For example, an accelerated rise in the trend component combined with an abnormally high amplitude in the periodic component may jointly indicate high risk. It can be understood that the synthesis of risk spectra is the process of compressing and refining complex, multi-dimensional spatiotemporal feature field information into key signal sequences that are directly related to geological hazard risk and are easier to understand and calculate.

[0044] Optionally, in the spatial interpolation stage, in addition to linear interpolation, Kriging interpolation based on geostatistics can be considered. Kriging interpolation can utilize the spatial variogram model of eigenvalues ​​to provide optimal unbiased estimates and estimates of interpolation errors. For time series decomposition, nonparametric methods such as singular spectrum analysis can be tried. Singular spectrum analysis does not assume a fixed period and can better capture complex time series structures. When identifying patterns, anomaly detection algorithms from machine learning, such as isolated forests or single-class support vector machines, can be introduced to automatically identify anomalous points or anomalous periods in the feature sequence that deviate from the normal pattern. The risk indicator set can be further integrated with external triggering factors, such as performing correlation analysis between the occurrence time of heavy rainfall events and feature mutation points to enhance the physical meaning of risk indicators. The data structure of the risk spectrum can be designed hierarchically to distinguish risk signals at different spatial scales (such as regional scale and local key point scale) to facilitate hierarchical risk assessment.

[0045] To comprehensively evaluate the contribution of different modes in a multi-resolution feature field to the overall risk, a method for measuring mode risk contribution can be defined. An exemplary formula for calculating mode risk contribution is as follows:

[0046] in: Representing the A risk contribution metric for each identified pattern (e.g., a mutation or a periodic oscillation), with a higher value indicating a greater potential risk contribution from the pattern. This represents the absolute change in the mean of the relevant feature sequences before and after the occurrence of the pattern, and is used to quantify the mutation intensity of the pattern. The standard deviation of the residual components of the feature sequence is used as a measure of the background noise level. This indicates the spatial extent affected by the pattern, which can be the total number of affected grid cells or their total area. Indicates the duration of the pattern. This represents a reference constant with the dimension of time. It is a small positive constant used to prevent the denominator from being zero when the pattern duration is extremely short, see Table 1.

[0047] Table 1: Risk Indicator Table ; Example 5: Configuring the parameters of the geological hazard risk calculation engine. The geological hazard risk calculation engine is a core calculation module based on numerical simulation methods. Its parameter configuration includes material property parameters and boundary condition parameters. The material property parameters are set for loess media. For example, the elastic modulus of loess ranges from 50 MPa to 200 MPa, the Poisson's ratio ranges from 0.25 to 0.35, the cohesion ranges from 20 kPa to 60 kPa, and the internal friction angle ranges from 25 degrees to 35 degrees. The boundary condition parameters include displacement constraints on the model region, such as setting the bottom of the model as a fixed constraint and the lateral side of the model as a normal roller constraint, as well as external load conditions, such as the self-weight stress field calculated based on the thickness of the overlying soil layer. The risk spectrum is input as an external excitation. The risk indicators contained in the risk spectrum are converted into excitation conditions acting on the geological body. For example, the "sudden drop in groundwater level" mode indicated in the risk spectrum is converted into a boundary condition for pore water pressure dissipation applied in the corresponding grid cell and time. The dissipation rate is proportional to the intensity value of the risk indicator. Numerical simulations are performed to calculate the stress, strain, or displacement response at each critical point of the geological hazard. The numerical simulations employ the finite element method or finite difference method, with the simulation time step consistent with the temporal resolution of the risk spectrum. At each time step, the system based on the coupling of equilibrium equations, constitutive relations, and geometric equations is solved, outputting the stress tensor, strain tensor, or displacement vector of the critical point of the geological hazard (e.g., the sliding surface of a potential landslide, the toe of a building foundation pit). Risk indicators are extracted from the response data. These risk indicators are quantitative values ​​derived from the numerical simulation results and used to directly measure stability. For example, a safety factor is calculated for each point based on stress results and material strength parameters; the safety factor is defined as the ratio of shear strength to shear stress. Alternatively, the maximum shear strain value is extracted as an indicator of the degree of deformation; or the magnitude and rate of change of the resultant displacement at the critical point are calculated. All risk indicators for the critical points are compiled to generate a risk state dataset. The risk state dataset is a structured table, with each row recording all risk indicator values ​​for a specific critical point of a geological hazard at a specific time point. Table fields include timestamp, critical point identifier, safety factor, maximum shear strain, displacement, and displacement rate.

[0048] Risk thresholds are set, which are critical values ​​determined based on engineering safety standards and historical disaster cases. For example, the risk threshold for the safety factor is set to 1.25, and the risk threshold for the displacement rate is set to 5 mm per day. The risk status dataset is traversed, marking risk values ​​exceeding the thresholds. The traversal process involves scanning the risk status dataset line by line. For each record, the value of each risk indicator is checked to see if it is greater than (or less than) its corresponding risk threshold for the safety factor. If any risk indicator exceeds the limit, an "exceeding limit" mark is added to the record, along with the specific indicator name and the value exceeding the limit. The number of times each key point is marked is counted as the event occurrence count. The statistical operation is grouped by geological disaster key points. For each key point, the total number of times it is marked as exceeding the limit across all time steps is calculated; this total count is the event occurrence count for that key point. The geological disaster risk knowledge base is queried. This knowledge base is a database storing historical disaster data, expert rules, and numerical simulation benchmark cases. A preset risk curve is retrieved. The risk curve is a function with the event occurrence count as the independent variable and the risk level as the dependent variable, typically existing as discrete data point pairs or piecewise functions. The process involves inputting the number of events occurring, retrieving the corresponding risk level from the risk curve, and for a specific geological hazard key point, using its calculated event occurrence count as input, interpolating or searching on the risk curve to obtain the corresponding risk level assignment. Risk levels are typically divided into a finite number of ordered categories such as low risk, medium risk, and high risk. A risk level is assigned to each key point, and the evaluation results are stored. The identifier of each geological hazard key point, the calculated event occurrence count, and the risk level assigned according to the risk curve are correlated to form the final risk assessment result table. This table can be output as a file or stored in a database, and can be further used to generate a risk level distribution map.

[0049] In some embodiments, the numerical simulation of the geological hazard risk calculation engine may need to consider the coupling effect of groundwater seepage and stress field, i.e., perform fluid-structure interaction analysis. In this case, the parameter configuration also needs to include hydrogeological parameters such as the permeability coefficient and specific yield of loess, and the governing equations also need to include the seepage equation and its coupling terms with the stress equation. The risk threshold may not be a single fixed value, but rather a dynamic threshold, for example, adjusting the displacement rate threshold according to seasonal precipitation patterns, and using a stricter threshold during the rainy season. The statistical method for the number of events can be refined, for example, distinguishing events that exceed the limits of different risk indicators and assigning weights to the severity of different indicators, thereby calculating the weighted number of events. The risk curves in the geological hazard risk knowledge base may not be globally uniform, but rather multiple risk curves are preset for different types of geological hazard key points or different geological conditions. When querying, it is necessary to select the applicable risk curve based on the attributes of the key point. It can be understood that the mapping from the physical response of the numerical simulation to the final risk level is a process of transforming continuous calculation results into discrete, decision-making-friendly categorical information, and the risk curve plays a crucial bridging role in this process.

[0050] Optionally, when performing numerical simulations, the initial stress field equilibrium can be considered. A stress equilibrium analysis can be performed first to ensure the model is in equilibrium in its initial state before applying dynamic excitation derived from the risk spectrum to improve simulation accuracy. For calculating risk indicators, in addition to the safety factor, the concept of reliability analysis can be introduced to calculate the failure probability as a risk indicator. This requires handling the spatial variability of material parameters. When counting the number of events, the concept of a time window can be introduced, counting only the number of events exceeding limits within the most recent period to reflect the current risk situation, rather than accumulating historical risks. The geological hazard risk knowledge base can be designed as a learning system, updating the risk curve to better reflect reality after a new hazard event occurs. The storage of risk assessment results can include confidence information to reflect the uncertainty of risk level assignment, for example, estimating the confidence level based on the slope of the risk curve within the interval of the event's occurrence frequency.

[0051] To quantify the relationship between the frequency of events and risk level in risk curve construction, a risk level transformation function can be defined. An example risk level transformation function is as follows:

[0052] in: Represented by the number of times the event occurs The resulting risk level code This indicates the number of times the event occurred, as statistically obtained. It is a small normal number used to smooth out the number of times the event occurs. The case where the logarithmic function is undefined when the logarithmic value is zero. It is a base parameter greater than 1, and its value determines the sensitivity of the risk level as the number of events increases. The smaller the value, the faster the level increases with the number of attempts.

[0053] See Figure 5 This chart is a core visualization result in the regional risk assessment conclusion generation process. It presents the number of key geological hazard points corresponding to different risk levels within the study area using a bar chart: the horizontal axis represents the risk level, and the vertical axis represents the number of key points, with colors corresponding to the risk levels. First, a risk calculation engine is used to obtain the risk indicators for each key point. Then, the dataset is scanned to count the number of events exceeding limits. Finally, the risk level is assigned by matching the risk curve from the risk knowledge base, and the number of key points at each level is counted. This chart intuitively reflects the overall distribution characteristics of regional risk. Low-risk key points account for the highest proportion, indicating that the overall regional risk is controllable. However, the existence of high-risk key points requires targeted and strengthened monitoring. It is a crucial visualization tool for integrating the risk levels of key points to form the regional assessment conclusion.

[0054] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

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

Claims

1. A method for assessing urban geological hazard risk in loess areas based on groundwater level, characterized in that, The method includes: A monitoring equipment network was deployed to continuously collect raw groundwater level data in urban areas of the Loess Plateau region; noise filtering and missing value imputation were performed on the raw groundwater level data to generate a clean monitoring dataset. By processing the clean monitoring dataset, an initial spatiotemporal sequence with a unified format is constructed through time alignment and spatial registration. Feature extraction techniques are then used to analyze the interaction pattern between groundwater level and geological parameters from the initial spatiotemporal sequence, and a coupled feature vector is output. The investigation examines the resolution distribution of the monitoring data and divides the region into sub-blocks of different levels based on the resolution. A customized grid specification, including cell size and grid density, is then established for each sub-block. Initiate the mesh generation process, create a multi-resolution mesh framework based on customized mesh specifications, and perform spatial allocation and temporal expansion of coupled feature vectors within the multi-resolution mesh framework to generate a multi-resolution feature field. Examine the time series of multi-resolution feature fields, extract risk indicators using pattern recognition methods, and combine them into a risk spectrum; The geological hazard risk calculation engine is invoked, and the risk spectrum is used as input to solve for the risk value of key geological hazard points, obtain the risk status dataset, scan the risk status dataset, detect risk exceeding the limit events of each key geological hazard point, and record the number of events. Access the geological disaster risk knowledge base, match risk curves based on the frequency of events, assign risk levels, integrate the risk levels of each key point, and form a regional risk assessment conclusion.

2. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The process of filtering noise and imputing missing values ​​from the raw groundwater level data to generate a clean monitoring dataset includes: The system receives raw groundwater level data streams, uses a sliding window algorithm to detect abnormal fluctuations, removes these abnormal fluctuations, and uses interpolation to fill in the data gaps. The processed data is standardized to adjust the data scale and distribution, making the data from all monitoring points comparable. Verify the integrity and consistency of the data, mark invalid data segments; remove invalid data segments, and finally output the clean monitoring dataset.

3. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The cleanroom monitoring dataset is processed by temporal alignment and spatial registration to construct an initial spatiotemporal sequence with a unified format, including: Read the cleanroom monitoring dataset, unify the time sampling frequency, and resample data from different frequencies to a common time axis; Align the spatial coordinate system to map the monitoring point locations to a standard geographic grid; extract the groundwater level time series for each grid point; Aggregate the time series of all grid points and organize them into an initial spatiotemporal sequence matrix according to time order and spatial location.

4. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The aforementioned feature extraction technique is used to analyze the interaction pattern between groundwater level and geological parameters from the initial spatiotemporal sequence, and output a coupled feature vector, including: Load the initial spatiotemporal sequence, combine it with the geological parameter database, including soil type and slope information; calculate the joint statistics of groundwater level and geological parameters; Principal component analysis is used to reduce the feature dimensionality while preserving the main direction of variation. Time sliding window analysis is then performed on the dimensionality-reduced features to enhance dynamic interaction signals. The enhanced feature set is then compiled to generate coupled feature vectors.

5. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The resolution distribution of the survey and monitoring data is used to divide the region into sub-blocks of different levels based on the resolution; a customized grid specification is provided for each sub-block, including cell size and grid density, including: Analyze the spatial resolution of the monitoring data and draw resolution contour maps; divide the region into high-resolution sub-blocks, medium-resolution sub-blocks, and low-resolution sub-blocks based on the contour maps; Set small cell size and high grid density for high-resolution sub-blocks; set medium cell size and medium grid density for medium-resolution sub-blocks; set large cell size and low grid density for low-resolution sub-blocks. Generate a grid specification configuration file for each sub-block.

6. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The initiated mesh generation process, which creates a multi-resolution mesh framework based on customized mesh specifications, includes: Parse the mesh specification configuration file and initialize the mesh generator; create mesh patches independently for each sub-block and record the patch vertices and connection relationships; Assemble the facets into a 3D mesh and assign volume element identifiers; check the continuity and boundary consistency of the mesh. Optimize the mesh structure, smooth the mesh transition areas, and output a multi-resolution mesh framework.

7. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, Within the multi-resolution mesh framework, the spatial allocation and temporal expansion of coupled feature vectors are performed to generate a multi-resolution feature field, including: Project the coupled feature vectors onto the vertices of the multi-resolution mesh frame; fill the feature values ​​within the mesh cells using a spatial interpolation algorithm; Extend the feature values ​​along the time axis to generate a feature snapshot at each time step; serialize all feature snapshots to construct a time series of a multi-resolution feature field.

8. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The time series of the multi-resolution feature field is examined, and risk indicators are extracted using pattern recognition methods and combined to form a risk spectrum, including: Decompose the time series of the multi-resolution feature field into trend components, periodic components, and residual components; calculate the statistical properties of each component; Identify patterns related to geological hazards, such as abrupt change points and periodic oscillations, quantify the intensity of these patterns, form a set of risk indicators, arrange the risk indicator set in chronological order, and synthesize a risk spectrum.

9. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The process involves calling the geological hazard risk calculation engine, taking the risk spectrum as input, solving for the risk values ​​of key geological hazard points, and obtaining a risk status dataset, including: Configure the parameters of the geological hazard risk calculation engine, including material properties and boundary conditions; input the risk spectrum as an external stimulus; Perform numerical simulations to calculate the stress, strain, or displacement response at each key point of the geological hazard; extract risk indicators from the response data. Compile risk indicators for all key points to generate a risk status dataset.

10. The method for assessing urban geological hazard risk in loess areas based on groundwater level as described in claim 1, characterized in that, The scanned risk status dataset detects risk exceeding events at each key geological hazard point and records the number of events occurring. Access the geological hazard risk knowledge base, match risk curves based on the frequency of event occurrences, and assign risk levels, including: Set a risk threshold, iterate through the risk status dataset, and mark risk values ​​that exceed the threshold; count the number of times each key point is marked as the number of events. Query the geological disaster risk knowledge base and retrieve the preset risk curve; input the number of times the event occurred, read the corresponding risk level from the risk curve, assign a risk level to each key point, and store the evaluation results.

Images