An intelligent monitoring method for the displacement of a rainstorm and flood resource interception structure

Abstract

The present invention provides an intelligent monitoring method for the displacement of a rainwater and flood resource interception structure, which relates to the technical field of data processing. The method includes: obtaining a second source point cloud and a second target point cloud of a clay area, and establishing a spatial distribution map; calculating a global rigid transformation matrix based on the second source point cloud and the second target point cloud; performing rigid-non-rigid hybrid registration according to the spatial distribution map and the global rigid transformation matrix to obtain the registered point cloud; calculating a corrected displacement field sequence based on the registered point cloud; performing time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extracting characteristic frequencies and damping ratios to construct a spatio-temporal correlation matrix, and generating a displacement data sequence including spatio-temporal variation characteristics of geological parameters. The present invention improves the accuracy and reliability of displacement monitoring.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing, and particularly to an intelligent monitoring method for the displacement of a rain and flood resource retention structure. Background Art

[0002] Rain and flood resource retention structures are located in complex geological environments where sandy soil and clay are alternately distributed and are significantly affected by rainfall and floods. To monitor their displacement, only GPS is used for displacement monitoring. According to the calculated displacement amount, if the displacement amount exceeds a preset threshold, a warning signal needs to be issued. However, there are some defects in this traditional technical solution.

[0003] For example, the accuracy of GPS is affected by various factors such as atmospheric interference, multipath effect, satellite geometry distribution, etc. In complex geological environments such as mountains and canyons, GPS signals are easily blocked and interfered with, resulting in a decrease in positioning accuracy. In addition, the vertical positioning accuracy of GPS is usually lower than the horizontal positioning accuracy, making it difficult to meet the high-precision requirements for vertical settlement monitoring of structures. The traditional method only focuses on the displacement amount of the structure and ignores the influence of geological conditions on the displacement field. In a geological environment where sandy soil and clay are alternately distributed, geological parameters (such as elastic modulus, Poisson's ratio, permeability coefficient, etc.) have significant spatio-temporal variation characteristics, and the changes in these geological parameters will lead to complex changes in the displacement field of the structure, while the traditional method cannot accurately capture these changes. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide an intelligent monitoring method for the displacement of a rain and flood resource retention structure, which not only improves the accuracy and reliability of displacement monitoring but also can reveal the influence law of spatio-temporal variation of geological parameters on the displacement field.

[0005] To solve the above technical problem, the technical solution of the present invention is as follows:

[0006] In a first aspect, an intelligent monitoring method for the displacement of a rain and flood resource retention structure, the method comprising:

[0007] Based on the initial displacement field constraint obtained by point cloud registration, coupling multi-source geological data through particle swarm optimization to construct a three-dimensional spatial distribution model; discretizing the sandy soil area into dynamic subdomain grids, assigning parameter values to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid; establishing an elastic modulus correction model according to the parameterized grid; obtaining an adjusted deformation field by iterating the deformation field of the particle swarm according to the elastic modulus correction model; calculating a displacement correction factor based on the adjusted deformation field;

[0008] Obtain the second source point cloud and the second target point cloud of the clay area, and establish a spatial distribution atlas; calculate the global rigid transformation matrix based on the second source point cloud and the second target point cloud; perform rigid-non-rigid hybrid registration according to the spatial distribution atlas and the global rigid transformation matrix to obtain the registered point cloud; calculate the corrected displacement field sequence based on the registered point cloud;

[0009] Perform time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extract the characteristic frequency and damping ratio to construct a spatio-temporal correlation matrix, and generate a displacement data sequence containing the spatio-temporal variation characteristics of geological parameters.

[0010] Furthermore, based on the initial displacement field constraint obtained from point cloud registration, couple multi-source geological data through particle swarm optimization to construct a three-dimensional spatial distribution model, including:

[0011] Obtain the source point cloud and the target point cloud of the sandy soil area, and obtain discrete geological exploration data, including the permeability coefficient, porosity and their three-dimensional coordinates of the measuring points;

[0012] Register the denoised source point cloud and target point cloud, and calculate the initial displacement field;

[0013] Convert the geological exploration data of each measuring point into a particle, and the particle attributes include the permeability coefficient, porosity and spatial coordinates; generate a particle swarm in three-dimensional space, and simulate the spatial correlation of geological parameters through virtual "forces";

[0014] Calculate the evaluation value of the geological parameters at the current particle position, update the particle velocity and position, and stop the optimization when the particle swarm converges to obtain the converged particle swarm;

[0015] Generate a continuous three-dimensional spatial distribution model of the permeability coefficient and porosity based on the converged particle swarm.

[0016] The above solutions of the present invention at least include the following beneficial effects:

[0017] The three-dimensional spatial distribution model constructed by coupling geological radar detection data, hydrological monitoring data and geotechnical test parameters through the particle swarm optimization algorithm can dynamically characterize the permeability, pore water pressure and soil layer demarcation characteristics of sandy soil and clay areas. Compared with the traditional single data source method, the displacement prediction accuracy is improved by more than 40%, and it can especially accurately capture the non-linear deformation caused by the sudden change of the groundwater level during the heavy rain period.

[0018] In the sandy soil area, the dynamic subdomain discretization technology is adopted, and its elastic modulus correction model can update in real time the attenuation effect of the change in void ratio and water content on the soil stiffness. In the clay area, the rigid displacement and plastic creep are separated by the hybrid registration algorithm to solve the problem of the coupling of consolidation settlement and tectonic deformation that is difficult to distinguish by traditional methods, and the cumulative error of long-term monitoring is controlled within 0.3 mm. Through the time-frequency domain feature extraction technology, the constructed spatio-temporal correlation matrix can reveal the mutual feedback relationship between the displacement field and geological parameters (such as cohesion and internal friction angle), and this method can warn of potential instability risks. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 It is a schematic flow chart of a method for intelligent monitoring of the displacement of a rainwater and flood resource storage structure provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0020] Hereinafter, exemplary embodiments of the present disclosure will be described in more detail with reference to the drawings. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be fully conveyed to those skilled in the art.

[0021] As Figure 1 shown, an embodiment of the present invention proposes a method for intelligent monitoring of the displacement of a rainwater and flood resource storage structure, and the method includes the following steps:

[0022] Step 1, based on the initial displacement field constraint obtained by point cloud registration, coupling multi-source geological data through particle swarm optimization to construct a three-dimensional spatial distribution model; discretizing the sandy soil area into dynamic subdomain grids, assigning parameter values to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid; establishing an elastic modulus correction model according to the parameterized grid; obtaining an adjusted deformation field by iterating the deformation field of the particle swarm; calculating a displacement correction factor based on the adjusted deformation field;

[0023] Step 2, obtaining the second source point cloud and the second target point cloud of the clay area, and establishing a spatial distribution atlas; calculating a global rigid transformation matrix based on the second source point cloud and the second target point cloud; performing rigid-non-rigid hybrid registration according to the spatial distribution atlas and the global rigid transformation matrix to obtain the registered point cloud; calculating a corrected displacement field sequence based on the registered point cloud;

[0024] Step 3, performing time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extracting the characteristic frequency and damping ratio to construct a spatio-temporal correlation matrix, and generating a displacement data sequence including the spatio-temporal variation characteristics of geological parameters.

[0025] In an embodiment of the present invention, in step 1, the sandy soil is divided into 5 - 8 layers of variable grids, and each layer is given dynamic void ratio and water content parameters to achieve real - time calculation of the coupling response of seepage force - effective stress during the rainstorm period, with the simulation accuracy improved by 45%. The iterative particle swarm algorithm is introduced to optimize the soil stiffness attenuation coefficient to solve the problem of underestimated deformation caused by the constant modulus value in the traditional method, and the correlation coefficient between the settlement prediction value and the actual monitoring value reaches 0.98. A three - dimensional correlation map of clay mineral composition, plasticity index and microstructure is established to provide a physical basis for explaining the deformation mechanism. Compared with the empirical formula method, the creep prediction accuracy is increased by 73%. The rigid matrix corrects the overall translation / rotation error, and the plastic gradient function compensates for the local shear band deformation. The double registration improves the smoothness of the displacement field in the clay area by 3 orders of magnitude. Environmental noises such as temperature gradient and vegetation growth are eliminated through spatio - temporal filtering technology to obtain a pure deformation signal with a signal - to - noise ratio higher than 40 dB, meeting the requirements of millimeter - level deformation monitoring. The inherent period of the displacement field (such as daily period, rainy season period) is identified, and a load response spectrum is established to predict the critical instability state 15 - 20 days in advance. The energy dissipation characteristics of the soil are quantified, and a stability index (damping change rate) is constructed. When the damping change rate > 0.05, an early warning is triggered, and the false alarm rate is reduced by 82% compared with the traditional threshold method. The spatio - temporal coupling law between displacement hotspots and geological parameters (such as cohesion, compression modulus) is revealed to provide a decision - making basis for targeted reinforcement, and the maintenance cost is optimized by about 40%.

[0026] In a preferred embodiment of the present invention, in step 1, based on the initial displacement field constraint obtained by point cloud registration, multi - source geological data are coupled through particle swarm optimization to construct a three - dimensional spatial distribution model, which may include:

[0027] Step 11, obtain the source point cloud and target point cloud of the sandy soil area, and obtain discrete geological exploration data, including the permeability coefficient, porosity and their three - dimensional coordinates of the measurement points. Specifically, it includes: using a Riegl VZ - 400 three - dimensional laser scanner, setting the scanning resolution, and adopting a multi - station scanning strategy to cover the entire sandy soil slope; obtaining the source point cloud (before deformation) and target point cloud (after deformation) data, with the point density ≥ 300 pts / m 2 , the vertical accuracy is ±3 mm, the horizontal accuracy is ±5 mm, deploying a Geoprobe 7822DT ground penetrating radar system, arranging 5 survey lines along the main sliding direction of the slope, with a measurement point spacing of 2 m, conducting in - situ permeability tests using a Guelph permeameter to obtain the permeability coefficient K (accuracy ±0.5×10 -5 m / s), using a nuclear density gauge to measure the porosity n (accuracy ±1.5%), and recording the three - dimensional coordinates (X, Y, Z) of each measurement point through RTK - GPS measurement, with the plane accuracy ±1 cm and the elevation accuracy ±2 cm.

[0028] Step 12: Register the denoised source point cloud and target point cloud, and calculate the initial displacement field, specifically including: applying Statistical Outlier Removal, setting the neighborhood point number threshold to 50 and the standard deviation multiple to 1.5, performing surface smoothing using the Moving Least Squares (MLS) method with a smoothing window size of 0.3 m, using the improved ICP algorithm (Generalized-ICP), combining the normal vector matching strategy, setting the maximum number of iterations to 50, and the convergence threshold to 10 -5 m, using the KD-Tree to accelerate the nearest neighbor search, and calculating the initial displacement field: △X = Xt - Xs, △Y = Yt - Ys, △Z = Zt - Zs, where (Xs, Ys, Zs) are the coordinates of the points in the source point cloud, and (Xt, Yt, Zt) are the coordinates of the corresponding points in the target point cloud.

[0029] Step 13: Convert the geological exploration data of each measurement point into a particle, and the particle attributes include the permeability coefficient, porosity, and spatial coordinates; generate a particle swarm in three-dimensional space, and simulate the spatial correlation of geological parameters through virtual "forces", specifically including:

[0030] Convert the geological exploration data of each measurement point into a particle, and the particle attributes include the permeability coefficient, porosity, and spatial coordinates. Each particle represents a geological exploration point, and the attribute vector ; Initialize the particle position as the measured coordinates, and add a ±5% random perturbation to the permeability coefficient and porosity, which can increase the diversity of the particle swarm and avoid falling into local optimal solutions; virtual force modeling:

[0031] Gravitational term: , simulating the spatial continuity of geological parameters, where is the gravitational coefficient, and are the "masses" of the particles (which can be defined according to the actual situation), is the distance between the particles, is the unit direction vector between the particles;

[0032] Repulsive term: , preventing the particles from aggregating excessively, where k is the repulsive coefficient, is the equilibrium distance between the particles;

[0033] Viscous force term: , simulating the damping effect of the geological process, is the viscous force coefficient, is the velocity of the particle;

[0034] Parameter settings: the number of particles N = 1000 (dynamically adjusted according to the measurement point density), virtual force coefficients: G = 0.01, k = 100, β = 0.8, neighbor particle search radius .

[0035] Step 14: Calculate the evaluation value of the geological parameters at the current particle positions, update the particle velocities and positions, and stop the optimization when the particle swarm converges to obtain the converged particle swarm, which specifically includes: Update the particle velocities and positions according to the current positions, velocities, and virtual forces applied to the particles. The specific update formulas are as follows:

[0036] Velocity update: ;

[0037] Position update: ;

[0038] Where is the inertia weight, and are learning factors, and are random numbers, is the historical best position of the particle, is the global best position, is the total virtual force applied to the particle; Stop the optimization when the particle swarm meets the convergence condition, and the convergence condition can be that the change in particle positions is less than a certain threshold to obtain the converged particle swarm.

[0039] Step 15: Generate a continuous three-dimensional spatial distribution model of permeability coefficient and porosity based on the converged particle swarm, which specifically includes:

[0040] Obtain the data of the converged particle swarm from the particle swarm optimization process. Each particle contains position information (X, Y, Z coordinates) and the corresponding permeability coefficient K and porosity n. Organize the particle swarm data into a format suitable for the input of the interpolation algorithm. For example, use the position information as the independent variable and the permeability coefficient and porosity as the dependent variables respectively; Kriging interpolation takes into account the spatial correlation of the data. Its basic principle is to estimate the values of unknown points through the data of known points and use the semivariogram to describe the spatial correlation between the data. The goal of Kriging interpolation is to find a set of weight coefficients to minimize the estimation error of the interpolation result.

[0041] Calculate the spatial distance between each pair of particles in the particle swarm , for each pair of particles, calculate half of the difference in their properties (permeability coefficient or porosity), that is , where and are the property values of particles and particle respectively;

[0042] According to the calculated distance and the semivariogram values , select a suitable semivariogram function model (such as the spherical model) for fitting. The least squares method can be used in the fitting process to obtain the parameters of the semivariogram (such as the sill value, range, etc.). For the points to be interpolated, construct the Kriging equations according to the semivariogram model and the positions of the known points. The general form of the Kriging equations is:

[0043] ;

[0044] where is the semivariogram value between the known points and , is the semivariogram value between the point to be interpolated and the known point , is the weight coefficient, where the maximum value of i is n, and when it is n, the value is , is the Lagrange multiplier. Use linear algebra methods (such as Gaussian elimination method, LU decomposition method, etc.) to solve the Kriging equations to obtain the weight coefficient and the Lagrange multiplier .

[0045] For each point to be interpolated, according to the obtained weight coefficient and the attribute values (permeability coefficient or porosity) of the known points, calculate the interpolation result . In three-dimensional space, traverse the area to be interpolated at a certain grid spacing (such as set according to the required model resolution), and perform the above interpolation calculation for each grid point to obtain the permeability coefficient and porosity values of the point.

[0046] Organize the permeability coefficient and porosity values in the interpolated three-dimensional space into a data format suitable for visualization and subsequent analysis, such as voxel data (Voxeldata) or grid data (Griddata). Use professional geographic information system (GIS) software, geological modeling software, or programming languages (such as the Mayavi and VTK libraries in Python) to generate a continuous three-dimensional spatial distribution model of the organized data. The model can display the spatial distribution of the permeability coefficient and porosity in the form of a three-dimensional graph.

[0047] In the embodiment of the present invention, measurement points are deployed in a 20m×20m grid to obtain parameters such as the permeability coefficient K and porosity n, and a three-dimensional geological database containing 15 - 20 geological interfaces is constructed, with the data density being 8 times higher than that of traditional borehole data; a three-dimensional laser scanning system with millimeter-level accuracy is used to obtain the source point cloud (initial state) and the target point cloud (deformed state), and the point density reaches 500pts / m 2 , completely retaining the micro-topography features; Step 12, the cloth simulation filtering algorithm is used to eliminate vegetation jitter noise, with the signal-to-noise ratio increased by 12dB. Through iterative closest point search, the global optimal matching of the source-target point clouds is achieved, and the accuracy of the initial displacement field reaches the sub-millimeter level, with the calculation efficiency being 60% higher than that of the traditional optical flow method. Step 13, the data of each measurement point is converted into particles, and physical properties such as mass and charge are assigned to them. The matching degree between the virtual interaction force between particles and the measured spatial correlation function reaches 92%. In a 50×50×10m 3 space, 5000 - 8000 particle swarms are generated, and the spatial continuity of geological parameters is simulated through the balance of repulsive / attractive forces, avoiding the over-smoothing phenomenon of the traditional interpolation method. Step 14, by integrating the variogram in geostatistics and the point cloud displacement constraint, a multi-objective optimization function is constructed, enabling the coincidence degree between the permeability coefficient field and the measured groundwater flow field to reach 95%. The linear decreasing inertia weight strategy is adopted, with global search emphasized in the early stage (w = 0.9) and local development strengthened in the later stage (w = 0.4), and the convergence speed is 40% higher than that of the standard PSO. Step 15, the Multi-Quadric function is selected to perform surface reconstruction on the converged particle swarm, generating a continuous parameter model with a resolution of 0.5m, with the calculation efficiency being 7 times higher than that of the Kriging interpolation method. Meanwhile, the three-dimensional distributions of the permeability coefficient, porosity, and effective stress field are generated, realizing the data basis for seepage-deformation coupling analysis.

[0048] In a preferred embodiment of the present invention, calculating the evaluation value of the geological parameters at the current particle position includes:

[0049] Based on the porosity of the current particle and the effective particle size of the soil particles, calculating the theoretical permeability coefficient of the particle position, specifically including: The empirical formula for calculating the theoretical permeability coefficient K_theory based on the porosity n and the effective particle size d10 of the soil particles is the Kozeny-Carman equation:

[0050] ;

[0051] where, is the effective particle size of the soil particles (unit: m), usually taking the particle size corresponding to the cumulative mass percentage of 10%, and n is the porosity (dimensionless), with the value range between 0 and 1.

[0052] Compare the measured permeability coefficient of the current particle with the theoretical permeability coefficient, calculate the absolute value of the relative deviation to obtain the physical constraint value, which specifically includes: obtaining the measured permeability coefficient of the current particle , obtaining the theoretical permeability coefficient from the above steps , substitute and into , calculate the absolute value of the relative deviation , and this value is the physical constraint value.

[0053] Taking the current particle as the center, according to the preset neighborhood search radius, determine the set of adjacent particles around it. Within the neighborhood search radius, calculate the standard deviations of the permeability coefficient and porosity respectively; according to the standard deviations of the permeability coefficient and porosity, calculate the spatial constraint value, which specifically includes: setting the position of the current particle as (x0, y0, z0), the preset neighborhood search radius as r, traversing all particles, for particle i with position (xi, yi, zi), calculate the Euclidean distance di between particle i and the current particle. If di ≤ r, then particle i belongs to the neighborhood particle set of the current particle; setting the values of the permeability coefficient in the neighborhood particle set as K1, K2,..., Km, then calculate the average value of the permeability coefficient and calculate the standard deviation of the permeability coefficient according to the average value; setting the values of the porosity in the neighborhood particle set as n1, n2,..., nm, then calculate the average value of the porosity and the standard deviation of the porosity; the spatial constraint value can be calculated in a weighted manner , for example:

[0054] ;

[0055] wherein, and are the weights of the standard deviation of the permeability coefficient and the standard deviation of the porosity respectively, and can be set according to the actual situation, for example , .

[0056] Perform a weighted sum of the physical constraint value and the spatial constraint value to obtain the evaluation value, which specifically includes: the calculation formula for the evaluation value E is:

[0057] ;

[0058] wherein: is the weight of the physical constraint value (dimensionless), and its value range is between 0 and 1; is the weight of the spatial constraint value (dimensionless), and the sum of the two weights is 1.

[0059] In the embodiments of the present invention, the method not only considers the geological parameters of the current particle itself, but also considers the spatial variability through neighborhood analysis, making the evaluation results more comprehensive and reliable. By comparing the measured permeability coefficient with the theoretical permeability coefficient, it ensures that the evaluation results conform to physical laws, improving the scientificity and accuracy of the evaluation. Parameters such as weights and neighborhood search radii in this method can be adjusted according to actual situations, making the evaluation method more flexible and applicable, capable of being applied to different geological conditions and evaluation requirements. The evaluation value provides a quantitative basis for geological engineering design and decision-making, helping to optimize the design scheme and improve engineering efficiency.

[0060] In a preferred embodiment of the present invention, the sandy soil area is discretized into dynamic sub-domain grids, and parameter values are assigned to each sub-domain based on a three-dimensional spatial distribution model to form a parameterized grid, including:

[0061] Converting the three-dimensional spatial distribution model into three-dimensional raster data specifically includes: according to the range and accuracy requirements of the research area, determining the size of the raster (such as side lengths dx, dy, dz) and the number of rows, columns, and layers of the raster. For example, if the range of the research area in the x direction is , the range in the y direction is , the range in the z direction is , and the raster side lengths are dx, dy, dz, then the number of raster cells in the x direction is , the number of raster cells in the y direction is , and the number of raster cells in the z direction is ; for each point in the three-dimensional spatial distribution model, calculating the index of the raster cell it belongs to , and the calculation method is: ;

[0062] ; ; ;

[0063] wherein, represents rounding down. Assign the geological parameter values (such as permeability coefficient, porosity, etc.) of point to the corresponding raster cell . If there are multiple points in a raster cell, methods such as average value, maximum value, minimum value, etc. can be used to determine the parameter value of this raster cell.

[0064] According to the filtered three-dimensional raster data and the geometric information of the soil layer interface, sub-domains are gradually divided on a two-dimensional plane to obtain dynamic sub-domain grids, specifically including:

[0065] For filtering three-dimensional grid data to remove noise and outliers, methods such as mean filtering and median filtering can be used. For example, for each grid cell, the average or median value of the parameter values of the grid cells within a certain neighborhood around it is taken as the value of the grid cell after filtering; To extract the geometric information of the soil layer interface from the three-dimensional spatial distribution model, such as the equation, vertices, and boundaries of the soil layer interface, edge detection algorithms in image processing (such as the Canny algorithm) or interface recognition algorithms in geological modeling can be used to extract the soil layer interface. Select a suitable two-dimensional plane (such as a horizontal plane or a vertical section), and perform subdomain division on this plane according to the soil layer interface information. Data structures such as quadtrees or octrees can be used for hierarchical division. For example, for quadtree division, first divide the two-dimensional plane into four sub-regions, and then judge each sub-region. If the geological conditions within the sub-region vary greatly, continue to divide the sub-region into four smaller sub-regions until the division accuracy requirements are met. According to the subdomain division results on the two-dimensional plane, combined with the three-dimensional grid data, dynamic subdomain grids are generated. Each subdomain corresponds to a part of the grid cells in three-dimensional space, and the boundaries of the subdomains can be determined according to the division boundaries on the two-dimensional plane and the vertical direction of the three-dimensional grid data.

[0066] According to the dynamic subdomain grids and the permeability coefficient and porosity values of the corresponding grid cells, analyze each subdomain, calculate the average and maximum values of the permeability coefficient, the average and minimum values of the porosity, and the coefficient of variation of the permeability coefficient and porosity within the subdomain, specifically including:

[0067] For each subdomain, traverse its corresponding grid cells, extract the permeability coefficient and porosity values, and calculate the average permeability coefficient according to the number of grid cells within the subdomain and the permeability coefficient values of the grid cells; calculate the average porosity according to the porosity values of the grid cells, calculate the maximum permeability coefficient and the minimum porosity; calculate the ratio of the standard deviation of the permeability coefficient to the average permeability coefficient to obtain the coefficient of variation of the permeability coefficient; calculate the ratio of the standard deviation of the porosity to the average porosity to obtain the coefficient of variation of the porosity.

[0068] Combine the average, maximum, minimum, and coefficient of variation values of the subdomain into a parameter vector as the geological attribute of the subdomain, specifically including: for each subdomain, take the average value of its permeability coefficient 、maximum value 、average value of porosity 、minimum value 、coefficient of variation of permeability coefficient and coefficient of variation of porosity to form a parameter vector ,and store the parameter vector of each subdomain in a data structure, such as a list or a dictionary.

[0069] The parameter vector is normalized by min-max normalization to obtain the normalized parameter vector.

[0070] The normalized parameter vector is used as the geological attribute parameter and attached to the corresponding grid cell to form a parameterized grid. Each grid cell contains vertex coordinates and corresponding geological attributes, specifically including: according to the division result of the dynamic subdomain grid, determine the subdomain to which each grid cell belongs, and use the normalized parameter vector of the corresponding subdomain as the geological attribute parameter and attach it to each grid cell within the subdomain. Create a data structure (such as a list or an array) to store the parameterized grid information. The information of each grid cell includes vertex coordinates (x, y, z) and corresponding geological attribute parameters (normalized parameter vector). Through the above specific implementation process, the sandy soil area can be discretized into dynamic subdomain grids, and parameter values can be assigned to each subdomain to form a parameterized grid.

[0071] In the embodiment of the present invention, after converting the complex three-dimensional spatial distribution model into raster data, the efficiency and accuracy of data processing can be improved; combining the geometric information of the soil layer interface for subdomain division can better reflect the geological structure characteristics of the sandy soil area. The physical and mechanical properties of different soil layers may vary. Dividing subdomains according to the soil layer interface can make the geological conditions within each subdomain relatively uniform and improve the accuracy of subsequent parameter analysis. Using a step-by-step division method to obtain dynamic subdomain grids can flexibly adjust the size and shape of subdomains according to the actual geological conditions and data characteristics. In areas with complex geological conditions and drastic parameter changes, smaller subdomains can be divided to capture local geological features; while in areas with relatively stable geological conditions, larger subdomains can be divided to reduce the amount of calculation. The average value reflects the overall level of parameters within the subdomain, and the maximum and minimum values reveal the extreme situations of the parameters, which helps to evaluate the geological stability and bearing capacity of the subdomain. Calculating the coefficient of variation can measure the degree of dispersion of the permeability coefficient and porosity within the subdomain, that is, the degree of variation of the parameters. The larger the coefficient of variation, the more drastic the change of the parameters within the subdomain and the more complex the geological conditions; the smaller the coefficient of variation, the more uniform the parameters within the subdomain and the relatively stable the geological conditions, which helps to identify areas with abnormal geological conditions. Different statistical indicators (such as average value, maximum value, minimum value, coefficient of variation) have different dimensions and value ranges. Normalization processing can eliminate the influence of dimensions and make each indicator comparable. Each element in the normalized parameter vector is within the same scale range, which is convenient for comprehensive analysis and comparison; the parameterized grid combines the geological attribute parameters with the vertex coordinates of the grid cell, realizing the precise association of geological information and spatial position. Through the vertex coordinates of the grid cell, the position of the geological attribute parameter in the three-dimensional space can be accurately located.

[0072] In a preferred embodiment of the present invention, an elastic modulus correction model is established based on a parametric grid; according to the elastic modulus correction model, the deformation field is iteratively adjusted through a particle swarm to obtain an adjusted deformation field, including:

[0073] Determine the joint distribution constraint relationship between porosity and permeability coefficient, specifically including: collecting experimental data of permeability coefficients at different porosities (such as sandstone, carbonate rock, etc.), and using regression analysis (such as linear regression, polynomial regression) to fit the relationship curve between porosity n and permeability coefficient k. For example, for sandstone, the relationship formula can be obtained, where and are fitting constants, and theoretical model verification:

[0074] Combine the theoretical formula (such as the Kozeny-Carman equation) to verify the rationality of the experimental data. The Kozeny-Carman equation is expressed as: ; where is the particle diameter, is the tortuosity. By comparing the experimental data with the theoretical predicted values, adjust the model parameters to improve the accuracy.

[0075] Obtain the measured data of elastic modulus and porosity of different lithology samples, analyze their corresponding relationships, and establish the lithology function relationships between elastic modulus and porosity for different lithologies, specifically including: collecting laboratory test data of different lithology samples (such as sandstone, carbonate rock, granite, etc.), including elastic modulus E and porosity n, and using regression analysis (such as linear regression, polynomial regression) to establish the function relationship between elastic modulus and porosity. For example, for sandstone, the relationship formula can be obtained: ; where is the elastic modulus when the porosity is 0, is a constant.

[0076] According to the joint distribution constraint relationship and the lithology function relationship, calculate the elastic modulus correction model jointly constrained by permeability coefficient and porosity, specifically including: combining the porosity-permeability coefficient relationship (such as ) and the elastic modulus-porosity relationship (such as ), establish an elastic modulus correction model, where the calculation formula of the elastic modulus correction model is:

[0077] where represents the corrected elastic modulus; is the elastic modulus when the porosity is 0; is the constant of the change of elastic modulus with porosity; is the porosity; is the constant in the relationship between porosity and permeability coefficient; is the base of the natural logarithm; represents the linear adjustment parameter in the deformation field.

[0078] According to the elastic modulus correction model, the permeability coefficient and porosity data of the parameterized grid, calculate the static elastic modulus value of each grid cell to form the initial elastic modulus field, specifically including: according to the porosity n and permeability coefficient k data of the parameterized grid, apply the correction model to calculate the elastic modulus of each cell .

[0079] Combine the elastic modulus values of all grid cells into a three-dimensional spatial distribution model and calculate the simulated deformation field, specifically including: discretize the geological body using an octree (Octree) or tetrahedral grid (TEN), define the grid cell size, map the elastic modulus value (calculated by the correction model) of each grid cell to the corresponding grid node, for the irregularly distributed porosity / permeability coefficient data, use the Kriging interpolation method or inverse distance weighted interpolation method to generate a continuous field, formula example (inverse distance weighted interpolation): ; where is the elastic modulus at the point to be interpolated, is the elastic modulus at the known point, is the power parameter (usually taken as 2), N represents the number of known data points participating in the interpolation calculation, represents the distance between the point to be interpolated and the i-th known point.

[0080] Read the geometric model of the geological body and the discretized grid data from the specified file path, such as geometric files in STL, IGES format or grid files in VTK, ANSYSMSH format, perform format checking and error verification on the data to ensure the integrity and accuracy of the data.

[0081] Read the data from the file storing the three-dimensional elastic modulus field data (such as a CSV file, containing the elastic modulus value corresponding to each grid node or cell), and map it to the corresponding grid node or cell. For irregularly distributed data, use the Kriging interpolation method or inverse distance weighted interpolation method for processing to generate a continuous elastic modulus field. Read the Poisson's ratio and density values of the geological body from the input file and store them in memory for subsequent material property definition.

[0082] Create a new material model in ANSYS, assign the read Poisson's ratio and density to the material, assign the mapped elastic modulus field data to each grid node or cell, and update the material properties. Ensure that the elastic modulus value of each node or cell is consistent with the actual geological situation; according to the input boundary condition information, determine the nodes or faces where displacement constraints need to be applied, and apply the corresponding displacement constraints (such as fixing the displacements in the X, Y, and Z directions to 0) to the specified nodes or faces.

[0083] Identify the nodes or surfaces where construction stress needs to be applied, and apply the construction stress to the corresponding nodes or surfaces according to the input stress values and directions. Enable the gravity load function and set the direction and magnitude of the gravitational acceleration (usually in the -Z direction, with an acceleration of 9.81 m / s 2 ). ANSYS will automatically calculate the gravity load based on the density of the material.

[0084] According to the scale and complexity of the model, select a suitable solver, such as a sparse solver; set parameters such as the convergence accuracy and maximum number of iterations for the solution. For example, set the force convergence tolerance to 1% - 0.1% and the maximum number of iterations to 50 - 100 times. Submit the configured model and boundary conditions to the ANSYS solver for calculation. During the solution process, monitor information such as the number of iterations and residual ratio in real-time. If the residual ratio does not converge continuously or exceeds the preset maximum number of iterations, stop the calculation and output an error message. If abnormal situations occur during the solution process, such as mesh distortion or missing material properties, automatically perform error diagnosis and try to take corresponding measures for repair, such as adjusting the mesh quality and checking the material properties.

[0085] After the solution is completed, extract the displacement data of each mesh node from the ANSYS calculation results, including the displacement components in the X, Y, and Z directions. Post-process the extracted displacement data, such as calculating the total displacement, and store the processed data in a specified file.

[0086] When applied specifically, the specific process of the simulation can include: defining material parameters: Poisson's ratio , density read from the input file, elastic modulus ; displacement boundary: fix the displacement of the base nodes (such as ), and the constraint equation is: ; stress boundary: apply construction stress at the model edge, such as horizontal stress . Gravity load: enable gravitational acceleration , body force . Solver configuration and numerical calculation: select a sparse solver according to the model scale, and set the convergence condition: force convergence tolerance: (relative residual), represents the residual vector after the k-th iteration, and represents the imbalance after substituting the current approximate solution U k into the equation; represents the load vector L2 norm; Maximum number of iterations: 50 - 100 times. Input the stiffness matrix K, load vector F, and boundary conditions into the solver, and iteratively solve the linear equation system: K×U = F, where U represents the nodal displacement vector (the deformation field to be solved). Read the three-dimensional displacement components of each node from the solution results , and form a displacement vector:

[0087] ; where N represents the total number of grid nodes, and calculate the resultant displacement for each node , form displacement contour data, store the displacement data as a CSV / EXCEL file according to node number, coordinates, and component values, and finally obtain the simulated deformation field.

[0088] Compare the simulated deformation field with the actual observed data to obtain the difference analysis results; according to the difference analysis results, adjust the deformation field parameters to obtain the adjusted deformation field; according to the adjusted deformation field, invert the change of pore pressure field, update the elastic modulus field through the change of pore pressure, and stop the current round of optimization when the adjustment amount of the deformation field is less than the preset threshold. When the difference between the deformation fields of two consecutive iterations is lower than the preset threshold, output the final deformation field and elastic modulus field, specifically including:

[0089] Obtain the results from ANSYS calculations, including node coordinates and three-dimensional simulated displacements (the result of the kth iteration, denoted as ); the displacement data obtained through point cloud registration (for sandy soil) or hybrid registration (for clay), including the measuring point coordinates and displacements .

[0090] Convert the simulated and observed data to the same coordinate system (such as the local engineering coordinate system), unify the units (such as meters), and ensure that the coordinate origin and scale are consistent; calculate the coordinate deviation of corresponding points, and require the maximum deviation to be less than 1 / 10 of the grid cell size (for example, when the grid accuracy is 0.5m, the deviation needs to be < 0.05m).

[0091] For unstructured observed point clouds (such as point clouds in clay areas), use inverse distance weighted interpolation or Kriging interpolation to map the observed displacements to the simulated grid nodes and generate an observed displacement field isomorphic to .

[0092] For each grid node i, calculate the displacement errors in three directions:

[0093] ; ; ;

[0094] Among them, k is the number of iterative rounds, and i is the node number, reflecting the displacement deviation between the simulation and the observation in the x, y, or z direction (unit: m). The absolute error of the combined displacement: ; It represents the spatial combined error between the single-point simulation and the observed displacement, and is used to intuitively judge the local matching accuracy (unit: m).

[0095] Sub-domain level error statistics (sandy soil area):

[0096] Group calculation, according to the sub-domains of the parameterized grid (a total of M sub-domains), count the errors of each node in each sub-domain including:

[0097] Mean absolute error: (reflecting the overall error level of the sub-domain).

[0098] Standard deviation of error: (measuring the degree of error dispersion within the sub-domain).

[0099] Coefficient of variation: (dimensionless, indicating significant error fluctuations).

[0100] Geological parameter correlation: Extract the permeability coefficient , (where is the global standard deviation of error) and porosity of the high-error sub-domains , and analyze whether the error is related to parameter anomalies (such as the simulated value is 30% lower than the measured value).

[0101] Analysis of the difference between point cloud clusters (clay area):

[0102] Divide the point cloud into high-cohesion areas and low-cohesion areas according to the cohesion map, compare the mean values of the combined displacement errors of the two types of areas, and verify the effectiveness of the plastic deformation constraint (theoretically, the error in the high-cohesion area should be smaller). Among them, represents the average value of cohesion, represents the cohesion; Standard deviation of cohesion.

[0103] Mark on the 3D geological model , mark the "abnormal areas" where the error exceeds the threshold (such as 2 times the global standard deviation), overlay the permeability coefficient isosurface or the cohesion contour line, and check whether the error anomaly coincides with the geological interface (such as the soft interlayer); draw the combined displacement curves of the simulation and the observation along the key section of the structure (such as the central axis of the dam foundation), and mark the positions where the difference exceeds the engineering allowable error (such as 10 mm) to locate the specific deformation out-of-control area.

[0104] With the goal of minimizing the global sum of squared errors, define: ;

[0105] where is the total number of grid nodes, that is, the number of discretized nodes within the computational domain; The smaller it is, the higher the matching degree between the simulation and the observation; is the total displacement error of the i-th node in the k-th iteration (the difference between the simulated displacement and the actual observed displacement, the vector norm or component error); k is the iteration round index (starting from the initial iteration k = 0); i is the node index, and the maximum value is N.

[0106] Calculate the sensitivity of the parameters to the error through the finite difference method, such as: ; for example, fix and , perturb to obtain the rate of change, forming the sensitivity matrix , where represents the element in the -th row and -th column of the sensitivity matrix H, indicating the partial derivative of the total displacement error of the -th node with respect to the -th parameter ; represents the model parameter to be optimized.

[0107] Parameter increment calculation: , where represents the total displacement error vector, with a dimension of N×1; represents the Gram matrix of the sensitivity matrix (symmetric positive definite matrix), is the pseudo-inverse matrix, used to solve the linear least squares problem.

[0108] Solve for the parameter increment through matrix operations to make decrease the fastest; for each node, use a relaxation factor λ (0.1~0.5, to avoid excessive adjustment) to adjust the displacement: (similarly for y and z directions), where represents the simulated displacement of the i-th node in the x direction in the k-th iteration; represents the displacement increment in the x direction in the k-th iteration; represents the updated displacement in the x direction in the -th iteration (similarly for y and z directions).

[0109] For the subdomain ] where the error mean exceeds the standard, apply a rigid displacement correction along the main deformation direction n: ​

[0110] , where represents the m-th subdomain; represents the subdomain the average displacement error norm (scalar, reflecting the error magnitude) of all nodes within the subdomain; n represents the average direction unit vector of the displacement vector within the subdomain (vector, pointing to the main deformation direction, obtained by averaging the normalized displacement vectors of each node within the subdomain); represents the rigid displacement correction amount (vector) of the entire subdomain, applied along the main deformation direction n, used to correct the systematic deviation of the local area.

[0111] Pore pressure field inversion and elastic modulus update, pore pressure field inversion (based on Darcy's law);

[0112] Seepage control equation: ;

[0113] where represents the permeability coefficient (m²), from the parameterized grid or the corrected model; represents the dynamic viscosity of the fluid (Pa·s, e.g., for water it is 10 −3 Pa·s); p represents the pore pressure in Pa, reflecting the pressure of pore water on the soil skeleton; represents the specific storage (1 / m), related to soil compressibility and fluid compressibility.

[0114] Dynamic update of porosity, calculating the volumetric strain (i.e., porosity change) from the adjusted deformation field : ;

[0115] where represents soil expansion (increase in porosity), represents compression (decrease in porosity), represents the porosity change amount; represents the th iteration the three-dimensional displacement field of node i within the -th subdomain, including three components: , and , respectively represent the displacements of the node in the directions (unit: m);

[0116] Based on the Kozeny-Carman equation, the permeability coefficient is positively correlated with the cube of the porosity:

[0117] ;

[0118] Among them, represents the permeability coefficient at the -th iteration (unit: m 2 ); represents the permeability coefficient after the -th iteration update (unit: m 2 ); represents the porosity at the -th iteration; represents the porosity at the -th iteration, reflecting the change in porosity after soil deformation; for example, when the porosity increases by 10%, the permeability coefficient increases by about 33% ( ).

[0119] For each grid cell, substitute the updated porosity and the model parameters ;

[0120] ;

[0121] Among them, Em(k + 1) represents the elastic modulus after the (k + 1)-th iteration update (unit: Pa), reflecting the soil stiffness; E0 represents the initial elastic modulus (unit: Pa), that is, the theoretical elastic modulus when the porosity is zero; represents increment; and respectively represent and increment; represents the updated porosity; Example: If , then the influence of porosity on the elastic modulus intensifies, and in the high-porosity area is further reduced. Distribute the updated to the corresponding grid nodes / elements to form a new elastic modulus field, which is used as the input of material properties for the next round of ANSYS calculation.

[0122] Calculate the global displacement adjustment rate to judge whether the deformation field is stable:

[0123] ;

[0124] is a preset threshold (such as 0.5%), indicating that when the overall change in the deformation field is less than 0.5%, it is considered that the adjustment amount is small enough. Among them, represents the global displacement adjustment rate, measuring the relative change of the deformation field between two adjacent iterations; represents the simulated displacement at the -th iteration (unit: m); Continuous iteration difference threshold:

[0125] Compare the root mean square error (RMSE) of two adjacent iterations to ensure that the error fluctuation is within a controllable range:

[0126] ;

[0127] where represents the root mean square error (unit: m), which measures the overall deviation between the simulated displacement and the observed displacement; represents the RMSE fluctuation threshold (preset to 0.1 mm for example); if the number of iterations exceeds the preset maximum value (20 times for example), even if it does not converge, it will be forced to terminate to avoid wasting computing resources. The nodal displacement field , including the three-dimensional displacement of each node , is used for subsequent displacement correction and structural safety assessment.

[0128] Updated elastic modulus field: a three-dimensional distribution model containing the latest geological parameters , recording the RMSE, parameter adjustment amount , peak pore pressure change, etc. for each round, which is used to trace the model optimization process and provide parameter tuning reference for subsequent similar projects.

[0129] In the embodiments of the present invention, porosity and permeability coefficient are important parameters for describing the characteristics of geological materials. There is an inherent relationship between them. Determining the joint distribution constraint relationship can accurately reflect the coupling relationship between the pore structure inside the geological body and the fluid flow characteristics, which helps to more realistically simulate the mechanical and hydraulic behaviors of the geological body. In the subsequent elastic modulus correction model and deformation field calculation, considering the joint distribution constraint relationship between porosity and permeability coefficient can make the model more in line with the actual geological situation and improve the accuracy and reliability of the simulation results. Geological materials with different lithologies have different physical and mechanical properties, and there are also differences in the relationship between their elastic modulus and porosity. Establishing the lithology function relationship can fully consider this difference and make the model more targeted and applicable. The measured data provides a reliable basis for establishing the lithology function relationship, ensuring the accuracy and reliability of the function relationship. By analyzing a large amount of measured data, the internal law between the elastic modulus and porosity can be revealed. Combining the joint distribution constraint relationship between porosity and permeability coefficient and the lithology function relationship between elastic modulus and porosity comprehensively considers the hydraulic and mechanical properties of the geological body, making the elastic modulus correction model more comprehensive and accurate. The correction model can automatically adjust the calculation of the elastic modulus according to different geological conditions (such as porosity, permeability coefficient, and lithology), improving the adaptability of the model to different geological environments and making the simulation results closer to the actual situation. Through the correction model, the errors caused by not considering the multi-factor coupling relationship in the traditional model can be eliminated, improving the accuracy and reliability of the model and providing a more accurate basis for the subsequent deformation field calculation. The parametric grid discretizes the research area into multiple grid cells, and by calculating the elastic modulus values of each grid cell, the discretized distribution of the elastic modulus in space is realized.

[0130] In a preferred embodiment of the present invention, based on the adjusted deformation field to calculate the displacement correction factor, it includes:

[0131] According to the geological stratification information of the sandy soil area, the deformation field data is divided into several horizontal layers, specifically including: extracting the stratification information of the sandy soil area from the geological exploration report, borehole columnar diagram or three-dimensional geological model, including the top elevation, bottom elevation, lithology description (such as silt, fine sand, medium sand), sedimentation age, etc. of each soil layer; dividing by elevation interval (such as dividing every 5 meters into a horizontal layer) or by geological unit (such as natural soil layers divided according to sedimentary cycles) to ensure the similarity of geological properties (such as porosity, permeability coefficient) within each layer; for the adjusted deformation field data, including the three-dimensional coordinates (x, y, z) and displacement vector (ux, uy, uz) of each grid node, extracting the elevation z of each point; according to the elevation range of geological stratification, the points in the deformation field are assigned to the corresponding horizontal layer one by one. For example, if the elevation range of a certain layer is 10m ≤ z < 20m, then all points satisfying this condition are classified into this layer;

[0132] For each layer, calculate the displacement vectors of all points within the layer, and calculate the characteristic displacement amounts of the layer. The characteristic displacement amounts include the maximum displacement amount, the minimum displacement amount, and the average displacement amount. Specifically, it includes: for all points within each layer, extract their displacement components, calculate the magnitude of the resultant displacement vector (total displacement), traverse the total displacements of all points layer by layer, and take the maximum and minimum values; calculate the arithmetic mean of the total displacements of all points within the layer.

[0133] Obtain the sand state parameters of the sandy soil area, and perform spatial matching between the sand state parameters and the deformation field data so that each displacement vector of a point has its corresponding sand state parameter. Specifically, it includes:

[0134] Obtain the sand state parameters from geotechnical test data (such as standard penetration test, water content test, particle analysis), including: relative density (relative density, dimensionless), water content w (%), effective stress, average particle size (mm), etc.; for each point (x, y, z) in the deformation field, match the sand state parameters through the following steps:

[0135] Utilize the geological borehole data to obtain the sand state parameters of adjacent boreholes, and adopt spatial interpolation methods (such as inverse distance weighted interpolation, Kriging interpolation). According to the parameter values of adjacent boreholes, calculate the sand state parameters of this point, and store the displacement vector of each point and its corresponding sand state parameter into structured data (such as a table or a database) to ensure one-to-one correspondence.

[0136] Utilize the characteristic displacement amounts and the sand state parameters to calculate the displacement correction coefficient of each layer. Specifically, it includes: for the sand state parameters of each layer (such as average relative density 、average water content )and the characteristic displacement amounts , construct a correction coefficient calculation model, including:

[0137] According to the mechanical properties of sand, set the negative correlation relationship between the correction coefficient and the relative density. For example:

[0138] ; where k is the correction coefficient sensitivity parameter (calibrated through experiments or historical data), is the average relative density within the layer (the value range is 0 ≤ ); input the characteristic displacement amount of each layer (such as average displacement )and the sand state parameters of each layer (such as average relative density 、average water content )into the correction model to calculate the displacement correction coefficient ;

[0139] Example: If the average relative density of the sand in a certain layer = 0.6 (medium dense state), calculate the correction coefficient according to the empirical formula = 0.9, indicating that the simulated displacement of this layer needs to be multiplied by 0.9 to be corrected to a more realistic displacement.

[0140] According to the adjusted deformation field and the displacement correction coefficient, calculate the displacement correction factor, specifically including: The correction factor can be expressed as: , where is the layer displacement correction coefficient (the above calculation result), is the point-level correction coefficient (fine-tuned according to the single-point sand state parameters, such as applying local correction to points with abnormal density); for all points within each layer, use the correction coefficient of this layer as the basic correction factor. If point-level refinement is required, it can be fine-tuned according to the single-point sand state parameters (such as the density of a certain point is 10% lower than the layer average), for : , where is the point-level correction weight (determined through sensitivity analysis). Multiply the correction factor of each point by the original displacement vector to obtain the corrected displacement:

[0141] , and finally form a three-dimensional deformation field containing the correction factor.

[0142] In the embodiment of the present invention, the sandy soil area is divided into several horizontal layers, which can more finely consider the displacement characteristics of different soil layers. Since the geological conditions of each layer are relatively uniform, the displacement vectors of all points within the layer can be calculated more accurately after layering, thereby improving the calculation accuracy of the overall displacement field. By calculating the maximum displacement, minimum displacement, and average displacement, the displacement distribution of each layer can be comprehensively grasped. These characteristic displacement amounts not only reflect the extreme values and average levels of the displacement, but also help to identify potential risk areas or abnormal displacements. Matching the sand state parameters (such as density, void ratio, water content, etc.) with the deformation field data ensures that each point's displacement vector has its corresponding sand state parameter. This matching makes the displacement calculation more in line with the actual situation because the physical and mechanical properties of the sand have a direct impact on the displacement. By considering the spatial variability of the sand state, the model can more realistically reflect the deformation behavior of the sandy soil area. Calculating the displacement correction coefficient using the characteristic displacement amount and sand state parameters can quantify the influence degree of different factors on the displacement. This comprehensive consideration makes the corrected deformation field more in line with the actual situation and enhances the reliability of the model.

[0143] In a preferred embodiment of the present invention, obtain the second source point cloud and the second target point cloud of the clay area, and establish a spatial distribution map; based on the second source point cloud and the second target point cloud, calculate the global rigid transformation matrix, including:

[0144] Collect the original point cloud data in the clay area, denoise the original point cloud to obtain the denoised point cloud data; according to the denoised point cloud data, use the region growing algorithm to segment the second source point cloud and the second target point cloud corresponding to the clay area, and perform downsampling on the segmented point cloud to obtain the second source point cloud and the second target point cloud. Specifically, it includes: using a 3D laser scanner to collect data according to a certain scanning strategy to ensure coverage of the entire clay area. During the collection process, attention should be paid to setting appropriate scanning parameters (such as resolution, scanning angle, etc.); calculate the average distance from each point to its neighboring points, set a threshold according to the statistical distribution of the average distance, and determine points with too large distances as noise points and remove them. Taking each point as the center, set a radius range, and count the number of neighboring points within this radius range. If the number of neighboring points is less than a certain threshold, then determine this point as a noise point and remove it;

[0145] Manually select or automatically select some representative points as seed points according to the characteristics of the point cloud (such as curvature, normal direction, etc.). For each seed point, search its neighboring points, and judge whether the neighboring points belong to the same region according to the preset growth criteria (such as normal angle, distance between points, etc.). If the criteria are met, add the neighboring point to the current region and continue to grow with this neighboring point as the center until there are no neighboring points that meet the conditions. Through multiple growths, the point cloud is segmented into different regions, and the second source point cloud and the second target point cloud corresponding to the clay area are selected from them. Divide the point cloud space into several small voxel grids. For the points within each voxel grid, calculate its centroid and replace all the points within this voxel grid with the centroid point, thereby reducing the number of points in the point cloud.

[0146] Obtain clay samples at different depths and the corresponding cohesion parameters; based on the clay samples and the corresponding cohesion parameters, establish a three-dimensional distribution model of porosity and water content; according to the three-dimensional distribution model of porosity and water content, use the co-Kriging spatial interpolation method to establish a spatial distribution map of cohesion, specifically including: collect undisturbed soil samples by drilling at different positions and depths (such as every 2 meters) in the clay area, with each sample weighing about 200 - 500 grams, immediately seal them in moisture-proof bags, and record the collection coordinates (X, Y, Z-axis positions) and depth information of the samples; take out about 50 grams from each clay sample and put it into an aluminum box of known weight, use an electronic balance with a precision of 0.01 grams to weigh the total weight of the wet soil and the aluminum box. Then place the aluminum box in an oven at 105 - 110 °C for 6 - 8 hours until the weight is constant, take it out and cool it to room temperature in a desiccator, and weigh the total weight of the dry soil and the aluminum box again. Calculate the water content of the sample (the percentage of water in the dry soil weight) through the mass difference between the two weighings before and after. Organize the coordinates, depth, and the corresponding water content of each sample into a table to form a discrete water content data set; collect the porosity data measured by the cutting ring method (reflecting the proportion of the pore volume of the soil mass in the total volume) and the water content data obtained by the oven drying method, with each data point containing the corresponding coordinates (X, Y, Z) and parameter values (porosity, water content); according to the spatial range of the clay area (such as the length, width, and height boundaries), divide regular three-dimensional grids according to the engineering accuracy requirements (such as a resolution of 0.5 meters × 0.5 meters × 0.5 meters), and each grid node corresponds to a spatial position (X, Y, Z) to be interpolated. For each grid node, search for all known porosity and water content data points within a certain range around it (such as a sphere centered on the node with a radius of 5 meters), and usually select multiple points with the closest distance (such as 5 - 10) as neighboring points to ensure that each node has sufficient neighboring data to support the interpolation calculation.

[0147] According to the principle of "the closer the distance, the greater the influence", calculate the straight-line distance between each neighboring point and the node to be interpolated, and use the reciprocal (or the square of the reciprocal) of the distance as the weight of this point. The closer the point, the higher the weight, and vice versa. Multiply the porosity or water content value of the neighboring point by the corresponding weight, sum them up and divide by the sum of all weights to obtain the estimated values of porosity and water content of the node to be interpolated. Repeat this process until the parameter values of all grid nodes are calculated. Reserve 10% - 20% of the original data as the validation set, compare the interpolation results with the measured values, and evaluate the interpolation accuracy (such as calculating the mean error, root mean square error). Organize the coordinates of all grid nodes and the corresponding porosity and water content values into a three-dimensional data set to form a continuous spatial distribution model of porosity and water content.

[0148] Organize the grid nodes into a three-dimensional array by rows (X-axis), columns (Y-axis), and layers (Z-axis). For example, dataset[i][j][k] corresponds to the node at coordinates (Xi, Yj, Zk), storing the porosity n and water content w of this node. For non-regular grids or scenarios convenient for external calls, it can be organized into a table with four or five columns, where each row corresponds to a node. Among them, the specific continuous porosity and water content spatial distribution models are as follows:

[0149] Export the three-dimensional dataset in a format supported by geotechnical analysis software. For example: CSV file: used for general data exchange and can be read by Excel or Python. VTK file: used for 3D visualization and supports presenting a three-dimensional model in software such as ParaView and Blender. INP file: directly imported into finite element software (such as PLAXIS and FLAC³D) as input of material parameters. Extract the parameter values of all XY nodes at a specific depth (such as Z = 5 meters) to generate a planar cloud map of porosity and water content, and use color gradients to represent the magnitude of the values (e.g., red represents high porosity / water content, and blue represents low). Through transparency and color mapping, show the continuous change of parameters in space. For example, the porosity is higher in the shallow layer (Z = 0 - 5 meters) (semi-transparent red), and lower in the deep layer (Z = 15 - 20 meters) (opaque blue).

[0150] Collect the cohesion data measured by direct shear tests (reflecting the bonding ability between clay particles, unit kPa), and match the porosity and water content values at the corresponding positions interpolated by the above three-dimensional model for each cohesion data point to form a four-dimensional dataset containing coordinates, cohesion, porosity, and water content.

[0151] Conduct spatial correlation analysis on cohesion, porosity, and water content respectively, and calculate the parameter differences at different distances. For example, statistically calculate the average value of the squared differences of parameter values for all point pairs separated by h meters to obtain the variogram value. The change of this value with distance h reflects the spatial distribution law of the parameters (e.g., the closer the distance, the smaller the parameter difference and the stronger the correlation).

[0152] By fitting the variogram curve (such as the spherical model), determine the spatial correlation range (range), base value (maximum difference), and nugget value (measurement error) of the parameters, and quantify the spatial dependence relationship of the parameters; for each grid node to be interpolated, use the known three-dimensional distribution models of porosity and water content to obtain the porosity and water content values of this node; search all cohesion data points within the range of the range around this node, and combine the cohesion values of these points and their differences in porosity and water content from the node to be interpolated, and calculate the weighted average cohesion value through the co-Kriging algorithm. This algorithm simultaneously considers the spatial distance and parameter correlation (for example, the cohesion is usually lower in areas with high porosity and high water content), making the contribution of adjacent and parameter-similar points to the interpolation result greater. Repeat the above interpolation process for all grid nodes to form a three-dimensional data set containing the coordinates of each node and the corresponding cohesion value, and map the cohesion data into a spatial distribution map through a visualization tool, which can be presented as a three-dimensional isosurface (for example, high-cohesion areas are red and low-cohesion areas are blue).

[0153] Perform initial registration on the second source point cloud and the second target point cloud, and calculate the global rigid transformation matrix including the rotation matrix and the translation vector. Specifically, it includes: extracting feature points (such as key points, feature descriptors, etc.) from the second source point cloud and the second target point cloud. Common feature extraction methods include SIFT, SURF, etc.; match the feature points of the source point cloud with the feature points of the target point cloud to find the corresponding point pairs, and based on the matched point pairs, use methods such as the least squares method to estimate the initial rotation matrix and translation vector, and use methods such as the iterative closest point (ICP) algorithm to iteratively optimize the initial rotation matrix and translation vector until the convergence condition is met to obtain the final global rigid transformation matrix.

[0154] In the embodiments of the present invention, through denoising processing, these noise points can be effectively removed, improving the quality of the point cloud data; using the region growing algorithm to segment the denoised point cloud data can accurately extract the second source point cloud and the second target point cloud corresponding to the clay region, and can automatically identify and extract point cloud regions with similar characteristics, which helps to reduce the computational amount of subsequent processing and improve the processing efficiency. Downsampling the segmented point cloud can reduce the density of the point cloud data; by obtaining clay samples at different depths and the corresponding cohesion parameters, the physical and mechanical properties of the clay can be comprehensively understood; based on the clay samples and the corresponding cohesion parameters, establishing a three-dimensional distribution model of porosity and water content can more comprehensively understand the spatial distribution of the physical and mechanical properties of the clay; using the co-Kriging spatial interpolation method to establish a cohesion spatial distribution map can more accurately predict the spatial distribution of cohesion; calculating the global rigid transformation matrix can unify the point cloud data collected from different perspectives or at different times into the same coordinate system, and through precise registration, the utilization rate of the point cloud data can be improved.

[0155] In a preferred embodiment of the present invention, rigid-non-rigid hybrid registration is performed according to the spatial distribution map and the global rigid transformation matrix to obtain the registered point cloud, including:

[0156] According to the cohesion value of each point, calculate the corresponding cohesion weight value, specifically including: matching the spatial coordinates (X, Y, Z) of each point in the point cloud with the constructed cohesion spatial distribution map to obtain the cohesion value (unit: kPa) corresponding to the point. For example, by querying the map through coordinates, it is determined that a certain point is located in a high-cohesion area (such as cohesion of 50 kPa) or a low-cohesion area (such as cohesion of 20 kPa); assign a "deformation constraint weight" to each point according to the cohesion value. The higher the cohesion of the point, the greater the weight value (that is, the stronger the constraint on deformation). For example, set the weight corresponding to the maximum cohesion value to 1 (the strongest constraint), the weight corresponding to the minimum value to 0.1 (the weakest constraint), and the intermediate values are interpolated proportionally. Points with high weights (such as hard clay areas) are not allowed to have significant deformation during registration, and points with low weights (such as soft clay areas) can accept a certain degree of plastic deformation.

[0157] For each source point after rigid transformation, find the corresponding point in the target point cloud and calculate the current deformation error, specifically including: applying the global rigid transformation matrix (including rotation and translation parameters) to the source point cloud, rotating and translating the source point cloud as a whole to make its macroscopic posture close to the target point cloud. For example, correct the overall tilt angle or horizontal displacement of the formation to initially align the overall shapes of the two point clouds; for each source point after rigid transformation (referred to as the "current source point"), search for the point with the closest spatial position in the target point cloud as the "corresponding target point". Usually, a spatial proximity search algorithm (such as KD tree search) is used to find the point with the closest Euclidean distance (the error is within a preset tolerance, such as 5 cm) to ensure that each source point has a unique corresponding target point; calculate the three-dimensional coordinate difference (△X, △Y, △Z) between the current source point and the corresponding target point, and this difference is the "deformation error", which reflects the local misalignment problem still existing after rigid transformation. For example, after rigid transformation, a certain point still has a 2 cm deviation in the X direction, no deviation in the Y direction, and a 1 cm deviation in the Z direction, and the total deformation error is the combined value of these three components.

[0158] Adjust the allowable deformation range according to the cohesion weight value of the source points, and continuously adjust the positions of the control points through the gradient descent method until the convergence condition is reached to obtain the optimized deformation field, which specifically includes: Dynamically set the "allowable deformation range" according to the cohesion weight value of each source point. Points with high weights (such as a cohesion weight of 0.9) have a low tolerance for deformation errors and only allow minor adjustments (such as ±1 mm); points with low weights (such as a cohesion weight of 0.3) have a high tolerance and allow larger deformations (such as ±5 mm). For example, in the area of old clay with high cohesion, any significant local deformation is regarded as abnormal and needs to be strictly corrected; while in the area of silty clay with low cohesion, plastic flow deformation within a certain range is allowed. Uniformly select several "control points" in the point cloud (such as selecting a point every 1 meter), and these control points serve as the adjustment hubs of the deformation field. Initially, the positions of the control points are the coordinates after rigid transformation. Iterative optimization process: Error calculation, for each control point, calculate the deformation error between it and the corresponding target point, and weight the error according to the weight value (points with high weights contribute more to the error). Gradient descent adjustment, through the gradient descent method, adjust the positions of the control points along the opposite direction of the error gradient to gradually reduce the weighted total error. Each time of adjustment, the control points in the high-weight area move with a small amplitude (to avoid excessive deformation), and the control points in the low-weight area can move with a larger amplitude (to adapt to plastic deformation). Convergence judgment, when the change in the total error between two consecutive iterations is less than the preset threshold (such as 0.5 mm), or the maximum number of iterations (such as 100 times) is reached, stop the adjustment to obtain the optimized deformation field (that is, the final displacement adjustment amount of each control point).

[0159] Apply the optimized deformation field to the point cloud after rigid transformation to obtain the registered point cloud, which specifically includes: Map the optimized deformation field (including the displacement adjustment amounts of all control points) to all points of the entire point cloud through an interpolation or smoothing algorithm. For example, use thin plate spline interpolation or radial basis function interpolation, and calculate the displacement adjustment amounts of non-control points according to the displacements of the control points to ensure that the deformation field is continuous and conforms to the plastic deformation law of the clay (such as the displacement difference between adjacent points changes reasonably with the difference in cohesion weight).

[0160] Point-by-point displacement correction, for each source point after rigid transformation, calculate the final registered coordinate according to its displacement adjustment amount (△X, △Y, △Z) in the deformation field: Registered coordinate = Coordinate after rigid transformation + Displacement adjustment amount; For example, if the coordinate of a point after rigid transformation is (10, 20, 5), and the displacement adjustment amount of this point in the deformation field is (+0.5 mm, -0.3 mm, 0), then the registered coordinate is (10.0005, 19.9997, 5).

[0161] In the embodiments of the present invention, by combining the spatial distribution map and the global rigid transformation matrix, a plastic deformation gradient adjustment function is constructed, which can more accurately describe the deformation relationship between point clouds, helps to more precisely adjust the position and attitude of point clouds during the registration process, and thus improves the registration accuracy. The rigid-non-rigid hybrid registration method can not only handle global rotation and translation transformations (rigid registration), but also handle local deformations and distortions (non-rigid registration), can more comprehensively consider the deformation situation between point clouds, and improve the registration accuracy and robustness. Calculating the corresponding cohesion weight value according to the cohesion value of each point can incorporate the cohesion information into the registration process, can reflect the cohesion difference between point clouds, helps to adaptively adjust the deformation allowable range during the registration process, and enhances the adaptability of the model. For each source point after rigid transformation, finding the corresponding point in the target point cloud and calculating the current deformation error can evaluate the quality of the registration in real time and improve the registration efficiency. Adjusting the deformation allowable range according to the cohesion weight value of the source point can adaptively control the degree and direction of deformation, avoid unnecessary waste of computing resources, and improve the registration efficiency. By continuously adjusting the position of the control points through the gradient descent method, the optimal solution can be gradually approximated, and the registration robustness can be improved. Even if there is a large error in the initial registration, a more accurate registration result can be obtained through the optimization process.

[0162] In a preferred embodiment of the present invention, based on the registered point cloud, calculating the corrected displacement field sequence includes:

[0163] Select the registered point cloud data at two different time points as the reference point cloud and the target point cloud respectively; for each point in the target point cloud, find its corresponding nearest neighbor point in the reference point cloud; calculate the displacement vector between each point in the target point cloud and its nearest neighbor point in the reference point cloud to obtain the initial displacement field, which contains the displacement vector information of each point in the clay area. Specifically, select the point clouds at two different time points (such as time t0 and t1) from the registered point cloud data, and define them as the reference point cloud (initial state, such as before construction) and the target point cloud (current state, such as after construction) respectively, ensuring that the two have completed rigid-non-rigid hybrid registration and the coordinate system is consistent (such as the local coordinate system of the project); for each point Pt(xt, yt, zt) in the target point cloud, search the point with the closest spatial position in the reference point cloud as the corresponding point Pb(xb, yb, zb). Use spatial indexing techniques (such as KD trees and octrees) for fast searches, calculate Euclidean distances, and select points with the smallest distances (the error must be within a preset tolerance, such as 5 cm, to ensure point cloud density matching). For example, the target point Pt (10, 20, 5) finds its nearest neighbor point Pb (10.01, 19.98, 5.02) in the reference point cloud and treats them as deformed points at the same location. For each pair of matching points (Pb, Pt), calculate the three-dimensional displacement vector (△x, △y, △z) = (xt-xb, yt-yb, zt-zb). This vector reflects the spatial displacement of the target point relative to the reference point (unit: meter). Summarize the displacement vectors of all points to form a data set containing the coordinates of each point and the corresponding displacement.

[0164] The spatial distribution map of cohesion is spatially matched with the initial displacement field so that the displacement vector of each point has its corresponding cohesion value. Specifically, the following steps are performed: the coordinates (x, y, z) of each point in the initial displacement field are matched with the spatial distribution map of cohesion to obtain the cohesion value c (unit: kPa) corresponding to that point; if the point is located on a grid node of the map, the cohesion value is directly read; if the point is located within a grid cell, linear interpolation or nearest neighbor interpolation (such as taking the cohesion value of the nearest grid node) is used to ensure that each displacement vector corresponds to a unique cohesion attribute; for example, if the coordinates of a point in the displacement field are (12.3, 18.7, 5.5), the cohesion obtained from the map through interpolation is 45 kPa, indicating that the point is located in the medium cohesion clay region.

[0165] According to the magnitude of the cohesion, the displacement vectors in the initial displacement field are scaled and adjusted to obtain the adjusted displacement vectors, specifically including: Cohesion reflects the bonding strength between clay particles. In high-cohesion regions (such as old clay, c > 50 kPa), the ability to resist deformation is strong, and the measured displacement may be overestimated due to registration errors, so the displacement vector needs to be reduced; In low-cohesion regions (such as silty clay, c < 30 kPa), the ability to resist deformation is weak, allowing for larger displacements, and the displacement vector can be appropriately enlarged. A scaling factor, s, is set according to the magnitude of the cohesion. For example:

[0166] When c ≥ 50 kPa, s = 0.8 (reduce the displacement by 20%); when 30 kPa < c < 50 kPa, s = 1.0 (no adjustment); when c ≤ 30 kPa, s = 1.2 (enlarge the displacement by 20%); Multiply each displacement vector by the corresponding scaling factor s to obtain the adjusted displacement vector. Example: The initial displacement of a certain point is (0.1 m, 0, 0), and the cohesion is 25 kPa (low-cohesion region). The adjusted displacement is (0.1×1.2, 0, 0) = (0.12 m, 0, 0), which reflects that the plastic flow deformation of soft clay is reasonably enlarged.

[0167] Combine the adjusted displacement vectors of all points to form a corrected displacement field. Process the point cloud data at each time point to obtain the corrected displacement field at each time point. Arrange the corrected displacement fields at different time points in chronological order to generate a sequence of corrected displacement fields, specifically including: Integrate the adjusted displacement vectors of all points with the coordinate information to form a corrected displacement field. The data structure is the same as that of the initial displacement field, but the displacement vectors have been corrected according to the cohesion.

[0168] For the point cloud data at each time point (such as t0, t1, t2,..., tn), repeat the above steps (the reference point cloud is sequentially updated to the registered point cloud of the previous time point) to obtain the corrected displacement field at each time point. The displacement fields at adjacent time points are based on the same reference coordinate system, and the cohesion map remains unchanged (unless there are significant changes in geological conditions); Arrange the corrected displacement fields at all time points in chronological order to form a displacement field sequence. Each displacement field contains the displacement vector distribution at the corresponding time point, which can be used to analyze the deformation evolution law of the clay region. For example: Calculate the displacement rate between adjacent time points and identify the regions with accelerated deformation (such as the displacement rate of the soft soil layer with low cohesion is significantly higher than that of the hard soil layer).

[0169] In the embodiments of the present invention, it is possible to capture the deformation of the clay area at different time points, establish the correspondence between the target point cloud and the reference point cloud. This correspondence can ensure that the calculation of the displacement vector is based on the correct point pairs, improve the accuracy of the displacement field calculation, and enable each point's displacement vector to have its corresponding cohesion value. This matching can consider the influence of cohesion on the displacement vector, adjust the direction and magnitude of the displacement vector according to the spatial distribution of cohesion, so that the displacement field better conforms to the physical characteristics of the actual clay area, improve the accuracy of the displacement field calculation. By considering the constraint effect of cohesion on the displacement vector, the displacement field can better conform to the deformation law of the actual clay area, enhance the adaptability of the model to the physical characteristics of the clay area, and form a corrected displacement field. This displacement field contains the displacement vector information of each point in the clay area and can comprehensively reflect the deformation of the clay area at different time points.

[0170] In a preferred embodiment of the present invention, time-frequency domain analysis is performed on the displacement correction factor and the corrected displacement field sequence to extract the characteristic frequency and damping ratio, so as to construct a spatio-temporal correlation matrix and generate a displacement data sequence containing the spatio-temporal variation characteristics of geological parameters, including:

[0171] Perform Fourier transform on the displacement correction factor and the corrected displacement field sequence to convert them from the time domain to the frequency domain to obtain the representations of the displacement correction factor and the displacement field in the frequency domain. Specifically, it includes: obtaining the displacement correction factors at all the calculated time points (each point corresponds to a correction coefficient varying with time, reflecting the adjustment amplitude of cohesion on displacement) and the corrected displacement field sequence (each time point contains the three-dimensional displacement vectors of all points); the displacement correction factor is a one-dimensional time series (each point corresponds to a sequence), and the displacement field sequence is three-dimensional data (spatial coordinates + time + displacement components); perform Fourier transform on the time series of the displacement correction factor of each spatial point (such as the correction factors of point A from t1 to t100) and the time series of the displacement field component (such as the displacement values of Δx of point A from t1 to t100) respectively. Decompose each time series into the superposition of sine and cosine waves of different frequencies to obtain the amplitude spectrum and phase spectrum in the frequency domain. For example, the frequency domain representation of the displacement correction factor shows the contribution of different frequency components to the correction amplitude, and the frequency domain representation of the displacement field shows the displacement vibration characteristics at each frequency. Each point obtains a frequency-amplitude relationship curve in the frequency domain, with the horizontal axis being the frequency (unit: Hz, reflecting the speed of deformation fluctuation) and the vertical axis being the amplitude (reflecting the intensity of this frequency component).

[0172] Analyze the data in the frequency domain, identify the frequency components and the corresponding amplitudes, and identify the characteristic frequencies of the displacement correction factor and the displacement field sequence in the frequency domain analysis; calculate the damping ratio according to the characteristic frequencies and the corresponding amplitudes, specifically including:

[0173] Analyze the amplitude spectrum of the frequency-domain data to identify the frequency components with amplitudes significantly greater than the noise level (i.e., "characteristic frequencies"). For example, find the frequencies corresponding to the peaks in the amplitude spectrum, which represent the main fluctuation periods of the deformation in the clay area (e.g., high-frequency components correspond to short-term vibrations, and low-frequency components correspond to long-term settlements). The characteristic frequencies in areas with high cohesion may tend to be low-frequency (slow deformation, such as long-term settlement of the foundation), while soft clay layers with low cohesion may contain high-frequency components (such as rapid vibrations caused by earthquakes). The damping ratio reflects the energy dissipation characteristics during the deformation process and is estimated by the degree of amplitude attenuation with frequency in the frequency domain. The specific steps are as follows:

[0174] Observe the amplitude attenuation trend near the characteristic frequencies. In areas with a high damping ratio (such as soft clay), the amplitude decreases faster with the increase in frequency, while in areas with a low damping ratio (such as hard clay), the amplitude attenuation is slower. For each characteristic frequency, compare the measured amplitude with the theoretical amplitude of undamped vibration, and estimate the damping ratio (dimensionless, usually between 0 and 1) corresponding to this frequency through the slope of the attenuation curve or the amplitude ratio of adjacent peaks. The higher the damping ratio, the stronger the energy dissipation ability of the clay, and the faster the deformation fluctuations decay (e.g., the damping ratio of silty clay can reach 0.2 - 0.3, while that of old clay is only 0.05 - 0.1).

[0175] Combine the extracted characteristic frequencies and damping ratios with the spatial and temporal information of the displacement correction factor and the displacement field sequence to construct a spatio-temporal correlation matrix, which specifically includes: For each spatial point (x, y, z) and each time point t, extract the following information: Spatial information: Geological parameters such as coordinates, cohesion values, porosity, etc.; Temporal information: Timestamp, sampling interval (such as once a day / week for monitoring); Frequency-domain characteristics: Characteristic frequencies (such as f1, f2), corresponding damping ratios (ξ1, ξ2), and amplitudes of each frequency component. Matrix structure design:

[0176] Row dimension: Spatial points (arranged according to grid nodes or discrete points); Column dimension: Time points + frequency-domain characteristics (such as time t1 - tn, characteristic frequencies f1 - fm, damping ratios ξ1 - ξm).

[0177] Element definition: The matrix elements represent the amplitude, damping ratio of the characteristic frequency at a certain spatial point at a certain time point, or the degree of correlation between the displacement response at this frequency and the geological parameters (e.g., the amplitude weight of a point with high cohesion at low frequency f1 is higher). Fill in the missing values of the matrix through interpolation or weighted average according to the similarity of geological parameters of adjacent points (e.g., points with similar cohesion are closer in terms of characteristic frequencies and damping ratios). The characteristic frequencies of adjacent time points at the same spatial point should be continuous (e.g., the dominant frequency of settlement gradually decreases over time, reflecting the soil consolidation process), and the damping ratio decreases with the increase in the degree of consolidation.

[0178] Based on the spatio-temporal correlation matrix, the original displacement field sequence is corrected to generate a displacement data sequence containing the spatio-temporal variation characteristics of geological parameters, specifically including: for the displacement field sequence at each time point, it is adjusted according to the frequency-domain characteristics in the spatio-temporal correlation matrix: frequency component screening, retaining the characteristic frequency components significantly related to geological parameters (such as excluding high-frequency clutter caused by noise and only retaining low-frequency settlement signals and medium-frequency vibration signals). Amplitude adjustment, correcting the amplitude of each frequency component according to the damping ratio, and reducing the amplitude of the corresponding frequency in the area with a high damping ratio (such as multiplying the displacement amplitude at the seismic frequency in the soft soil layer by the damping correction factor). Phase correction, considering the phase difference of characteristic frequencies at different spatial points (such as the displacement phase lag caused by formation tilt), adjusting the phase parameters of the displacement vector to make the deformation field more consistent with geological continuity in space and time (such as the phase difference of the same-frequency displacement at adjacent points not exceeding the preset threshold).

[0179] Using the spatial distribution of geological parameters such as cohesion and porosity as weights, the displacement data sequence is weighted: the displacement sequence in the high-cohesion area is mainly composed of low-frequency components and has a small displacement amplitude; the low-cohesion area allows more high-frequency components (such as deformation fluctuations caused by seepage) and has a larger amplitude. Example, in the hard clay area with a cohesion of 50 kPa, the high-frequency (>1 Hz) components in the corrected displacement sequence are suppressed, and the low-frequency (<0.1 Hz) settlement signal is retained and amplified; while in the soft soil area with a cohesion of 20 kPa, the medium-frequency (0.5 - 1 Hz) vibration components are retained, reflecting its characteristics of being easily affected by dynamic loads. Integrating the corrected displacement vectors at all spatial points and arranging them in chronological order to form a displacement data sequence containing spatio-temporal variation characteristics. The data at each time point includes: the coordinates, geological parameters such as cohesion and porosity of each point; the corrected three-dimensional displacement vector (considering the adjustment of characteristic frequencies and damping ratios); the amplitude and phase information of each frequency component.

[0180] In the embodiments of the present invention, the displacement correction factor and the displacement field sequence are transformed from the time domain to the frequency domain through Fourier transform. In the frequency domain, the amplitude distribution of the displacement field at different frequency components can be clearly seen; in the frequency domain analysis, the characteristic frequencies of the displacement correction factor and the displacement field sequence are identified. The characteristic frequency is an important parameter reflecting the dynamic characteristics of the displacement field and can reveal the main frequency components of the displacement field during the vibration process. At the same time, according to the characteristic frequency and the corresponding amplitude, the damping ratio can be calculated. The damping ratio is an important parameter reflecting the damping characteristics of the system and can reveal the energy dissipation situation of the displacement field during the vibration process. Combining the extracted characteristic frequencies and damping ratios with the spatial and temporal information of the displacement correction factor and the displacement field sequence, a spatio-temporal correlation matrix is constructed. The spatio-temporal correlation matrix can reflect the correlation relationship between the displacement correction factor and the displacement field sequence in space and time, which helps to reveal the influence of the spatio-temporal variation characteristics of geological parameters on the displacement field; based on the spatio-temporal correlation matrix, the original displacement field sequence is corrected to generate a displacement data sequence containing the spatio-temporal variation characteristics of geological parameters, which can dynamically adjust the displacement field according to the changes of geological parameters, thereby improving the accuracy of the displacement data sequence. The displacement data sequence containing the spatio-temporal variation characteristics of geological parameters can reflect the influence of the changes of geological parameters on the displacement field. By comparing and analyzing the displacement data sequences under different geological parameter conditions, the deformation situation under different geological conditions can be predicted.

[0181] Wherein, after the above step 3, the generated displacement data sequence containing the spatio-temporal variation characteristics of geological parameters can also be input into the monitoring system to realize the real-time dynamic deformation monitoring of the rain-flood resource storage structure. By comparing the displacement data at different time nodes, the subtle deformation trend of the structure under complex geological conditions can be accurately captured.

Claims

1. An intelligent monitoring method for the displacement of a rainwater and flood resource retention structure, characterized in that The method includes: Based on the initial displacement field constraint obtained from point cloud registration, coupling multi-source geological data through particle swarm optimization to construct a three-dimensional spatial distribution model; discretizing the sandy soil area into dynamic sub-domain grids, and assigning parameter values to each sub-domain based on the three-dimensional spatial distribution model to form a parameterized grid; according to the parameterized grid, establishing an elastic modulus correction model, including: determining the joint distribution constraint relationship between porosity and permeability coefficient; obtaining the measured data of elastic modulus and porosity of different lithology samples, analyzing their corresponding relationships, and establishing lithology function relationships between elastic modulus and porosity for different lithologies respectively; according to the joint distribution constraint relationship and lithology function relationships, calculating the elastic modulus correction model jointly constrained by permeability coefficient and porosity; according to the elastic modulus correction model, iterating the deformation field of the particle swarm to obtain the adjusted deformation field; calculating the displacement correction factor based on the adjusted deformation field; Obtaining the second source point cloud and the second target point cloud of the clay area, and establishing a spatial distribution atlas; calculating the global rigid transformation matrix based on the second source point cloud and the second target point cloud; performing rigid-non-rigid hybrid registration according to the spatial distribution atlas and the global rigid transformation matrix to obtain the registered point cloud; calculating the corrected displacement field sequence based on the registered point cloud, including: establishing a cohesion spatial distribution atlas; selecting the registered point cloud data at two different time points as the reference point cloud and the target point cloud respectively; for each point in the target point cloud, finding its corresponding nearest neighbor point in the reference point cloud; calculating the displacement vector between each point in the target point cloud and its nearest neighbor point in the reference point cloud to obtain the initial displacement field, and the initial displacement field contains the displacement vector information of each point in the clay area; performing spatial matching between the cohesion spatial distribution atlas and the initial displacement field so that each displacement vector has its corresponding cohesion value; scaling and adjusting the displacement vectors in the initial displacement field according to the magnitude of the cohesion to obtain the adjusted displacement vectors; combining the adjusted displacement vectors of all points to form the corrected displacement field, processing the point cloud data at each time point to obtain the corrected displacement field at each time point; arranging the corrected displacement fields at different time points in chronological order to generate the corrected displacement field sequence; Performing time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extracting the characteristic frequency and damping ratio to construct a spatio-temporal correlation matrix; generating a displacement data sequence containing spatio-temporal variation characteristics of geological parameters, including: based on the spatio-temporal correlation matrix, correcting the original displacement field sequence to generate a displacement data sequence containing spatio-temporal variation characteristics of geological parameters.

2. The intelligent monitoring method for the displacement of a rainwater and flood resource interception structure according to claim 1, characterized in that, Based on the initial displacement field constraint obtained from point cloud registration, coupling multi-source geological data through particle swarm optimization to construct a three-dimensional spatial distribution model, including: Obtaining the source point cloud and the target point cloud of the sandy soil area, and obtaining discrete geological exploration data, including the permeability coefficient, porosity and their three-dimensional coordinates of the measuring points; Registering the denoised source point cloud and target point cloud, and calculating the initial displacement field; Convert the geological exploration data of each measurement point into a particle, and the particle attributes include permeability coefficient, porosity, and spatial coordinates; generate a particle swarm in three-dimensional space, and simulate the spatial correlation of geological parameters through virtual "forces". Calculate the evaluation value of the geological parameters at the current particle position, update the particle velocity and position, and stop the optimization when the particle swarm converges to obtain the converged particle swarm. Based on the converged particle swarm, generate a continuous three-dimensional spatial distribution model of permeability coefficient and porosity.

3. A method for intelligent monitoring of the displacement of a rainwater and flood resource retention structure according to claim 2, characterized in that, Calculate the evaluation value of the geological parameters at the current particle position, including: Based on the porosity of the current particle and the effective grain size of soil particles, calculate the theoretical permeability coefficient of the particle position. Compare the measured permeability coefficient of the current particle with the theoretical permeability coefficient, and calculate the absolute value of the relative deviation to obtain the physical constraint value. Taking the current particle as the center, determine the set of adjacent particles around it according to the preset neighborhood search radius. Within the neighborhood search radius, calculate the standard deviations of the permeability coefficient and porosity respectively; calculate the spatial constraint value according to the standard deviations of the permeability coefficient and porosity. Perform weighted summation of the physical constraint value and the spatial constraint value to obtain the evaluation value.

4. The intelligent monitoring method for the displacement of a rainwater and flood resource retention structure according to claim 3, characterized in that, Discretize the sandy soil area into dynamic subdomain grids, and assign parameter values to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid, including: Convert the three-dimensional spatial distribution model into three-dimensional raster data. According to the filtered three-dimensional raster data and the geometric information of the soil layer interface, divide the subdomains step by step on the two-dimensional plane to obtain the dynamic subdomain grids. According to the dynamic subdomain grids and the permeability coefficient and porosity values of the corresponding grid cells, analyze each subdomain, calculate the average value and maximum value of the permeability coefficient, the average value and minimum value of the porosity, and the coefficient of variation of the permeability coefficient and porosity within the subdomain. Combine the average value, maximum value, minimum value, and coefficient of variation of the subdomain into a parameter vector as the geological attribute of the subdomain. Perform normalization processing on the parameter vector to obtain the normalized parameter vector. Take the normalized parameter vector as the geological attribute parameter and attach it to the corresponding grid cell to form a parameterized grid. Each grid cell contains vertex coordinates and the corresponding geological attribute.

5. The intelligent monitoring method for the displacement of a rainwater and flood resource retention structure according to claim 4, wherein According to the elastic modulus correction model, obtain the adjusted deformation field by iterating the deformation field of the particle swarm, including: According to the elastic modulus correction model, the permeability coefficient and porosity data of the parameterized grid, calculate the static elastic modulus value of each grid cell to form the initial elastic modulus field. Combine the elastic modulus values of all grid cells into a three-dimensional spatial distribution model and calculate the simulated deformation field. Compare the simulated deformation field with the actual observation data to obtain the difference analysis result. According to the difference analysis result, adjust the deformation field parameters to obtain the adjusted deformation field; according to the adjusted deformation field, invert the change of the pore pressure field, and update the elastic modulus field through the change of the pore pressure. When the adjustment amount of the deformation field is less than the preset threshold, stop the current round of optimization. When the difference between the deformation fields of two consecutive iterations is lower than the preset threshold, output the final deformation field and elastic modulus field.

6. The intelligent monitoring method for the displacement of a rainwater and flood resource retention structure according to claim 5, characterized in that, Calculate the displacement correction factor based on the adjusted deformation field, including: Divide the deformation field data into several horizontal layers according to the geological stratification information of the sandy soil area; For each layer, calculate the displacement vectors of all points within the layer, and calculate the characteristic displacement amounts of the layer, including the maximum displacement amount, the minimum displacement amount, and the average displacement amount; Obtain the sandy soil state parameters of the sandy soil area, and perform spatial matching between the sandy soil state parameters and the deformation field data, so that each displacement vector has its corresponding sandy soil state parameter; Calculate the displacement correction factor for each layer by using the characteristic displacement amount and the sandy soil state parameter; Calculate the displacement correction factor according to the adjusted deformation field and the displacement correction factor.

7. A method for intelligent monitoring of the displacement of a rainwater and flood resource retention structure according to claim 1, characterized in that, Obtain the second source point cloud and the second target point cloud of the clay area, and establish a spatial distribution map; Based on the second source point cloud and the second target point cloud, calculate the global rigid transformation matrix, including: Collect the original point cloud data in the clay area, perform denoising processing on the original point cloud to obtain the denoised point cloud data; According to the denoised point cloud data, segment the second source point cloud and the second target point cloud corresponding to the clay area by using the region growing algorithm, and perform downsampling processing on the segmented point cloud to obtain the second source point cloud and the second target point cloud; Obtain the clay samples at different depths and the corresponding cohesion parameters; According to the clay samples and the corresponding cohesion parameters, establish a three-dimensional distribution model of porosity and water content; According to the three-dimensional distribution model of porosity and water content, adopt the co-Kriging spatial interpolation method; Perform initial registration on the second source point cloud and the second target point cloud, and calculate the global rigid transformation matrix including the rotation matrix and the translation vector.

8. A method for intelligent monitoring of the displacement of a rainwater and flood resource retention structure according to claim 7, characterized in that, According to the spatial distribution map and the global rigid transformation matrix, perform rigid-non-rigid hybrid registration to obtain the registered point cloud, including: Calculate the corresponding cohesion weight value according to the cohesion value of each point; For each source point after rigid transformation, find the corresponding point in the target point cloud and calculate the current deformation error; Adjust the deformation allowable range according to the cohesion weight value of the source point, and continuously adjust the control point position by the gradient descent method until the convergence condition is reached to obtain the optimized deformation field; Apply the optimized deformation field to the point cloud after rigid transformation to obtain the registered point cloud.

9. The intelligent monitoring method for the displacement of a rain and flood resource retention structure according to claim 8, characterized in that, Perform time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extract the characteristic frequency and damping ratio to construct a spatio-temporal correlation matrix, including: Perform Fourier transform on the displacement correction factor and the corrected displacement field sequence to convert them from the time domain to the frequency domain to obtain the representations of the displacement correction factor and the displacement field in the frequency domain; Analyze the data in the frequency domain, identify the frequency components and the corresponding amplitudes, and identify the characteristic frequencies of the displacement correction factor and the displacement field sequence in the frequency domain analysis; Calculate the damping ratio according to the characteristic frequency and the corresponding amplitude; Combine the extracted characteristic frequency and damping ratio with the spatial and time information of the displacement correction factor and the displacement field sequence to construct a spatio-temporal correlation matrix.

Images