Intelligent monitoring method for displacement of rainfall flood resource retaining structure

Abstract

The invention provides a rainfall flood resource retaining structure displacement intelligent monitoring method, and relates to the technical field of data processing, and the method comprises the steps: obtaining a second source point cloud and a second target point cloud of a clay region, and building 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 mixed registration according to the spatial distribution map and the global rigid transformation matrix to obtain a registered point cloud; calculating a corrected displacement field sequence based on the registered point cloud; and carrying out time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extracting a characteristic frequency and a damping ratio to construct a space-time incidence matrix, and generating a displacement data sequence containing geological parameter space-time variation characteristics. According to the invention, the precision and reliability of displacement monitoring are improved.

Description

Technical Field

[0001] The invention relates to the technical field of data processing, and in particular to an intelligent monitoring method for displacement of a rainwater and flood resource interception and storage structure. Background Art

[0002] The stormwater storage structures are located in a complex geological environment with alternating distribution of sandy soil and clay, and are significantly affected by rainfall and floods. In order to monitor their displacement, only GPS is used for displacement monitoring. According to the calculated displacement, if the displacement exceeds the preset threshold, an early warning signal needs to be issued. However, this traditional technical solution has some defects.

[0003] For example, the accuracy of GPS is affected by many factors, such as atmospheric interference, multipath effects, and satellite geometric distribution. In complex geological environments, such as mountainous areas and canyons, GPS signals are easily blocked and interfered with, resulting in reduced positioning accuracy. In addition, the vertical positioning accuracy of GPS is usually lower than the horizontal positioning accuracy, which makes it difficult to meet the high-precision requirements for vertical settlement monitoring of structures. Traditional methods only focus on the displacement of structures, while ignoring the impact 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 spatiotemporal variation characteristics. Changes in these geological parameters will lead to complex changes in the displacement field of the structure, and traditional methods 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 stormwater resource interception and storage structure, which not only improves the accuracy and reliability of displacement monitoring, but also can reveal the influence of the temporal and spatial variation of geological parameters on the displacement field.

[0005] In order to solve the above technical problems, the technical solution of the present invention is as follows: In a first aspect, a method for intelligently monitoring the displacement of a rainwater resource interception and storage structure is provided, the method comprising: Based on the initial displacement field constraints obtained by point cloud registration, a three-dimensional spatial distribution model is constructed by coupling multi-source geological data through particle swarm optimization; the sandy soil area is discretized into dynamic subdomain grids, and parameter values ​​are assigned to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid; based on the parameterized grid, an elastic modulus correction model is established; based on the elastic modulus correction model, the deformation field of the particle swarm is iterated to obtain the adjusted deformation field; based on the adjusted deformation field, the displacement correction factor is calculated; Obtain the second source point cloud and the second target point cloud of the clay area and establish a spatial distribution map; 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 map and the global rigid transformation matrix to obtain the registered point cloud; calculate the corrected displacement field sequence based on the registered point cloud; The displacement correction factors and the corrected displacement field sequences are analyzed in the time-frequency domain to extract characteristic frequencies and damping ratios to construct a spatiotemporal correlation matrix and generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters.

[0006] Furthermore, based on the initial displacement field constraints obtained by point cloud registration, a three-dimensional spatial distribution model is constructed by coupling multi-source geological data through particle swarm optimization, including: Obtain source point cloud and target point cloud of sandy soil area, and obtain discrete geological exploration data, including permeability coefficient, porosity and 3D coordinates of the measuring points; Align the denoised source point cloud and target point cloud and calculate the initial displacement field; The geological exploration data of each measuring point is converted into a particle, and the particle attributes include permeability, porosity and spatial coordinates; a particle group is generated in three-dimensional space, and the spatial correlation of geological parameters is simulated through virtual "force"; Calculate the evaluation value of the geological parameters of the current particle position, update the particle speed and position, and stop the optimization when the particle group converges to obtain the converged particle group; Based on the converged particle swarm, a continuous three-dimensional spatial distribution model of permeability and porosity is generated.

[0007] The above solution of the present invention includes at least the following beneficial effects: By coupling geological radar detection data, hydrological monitoring data and geotechnical test parameters with the particle swarm optimization algorithm, the constructed three-dimensional spatial distribution model can dynamically characterize the permeability, pore water pressure and soil boundary 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%, especially the nonlinear deformation caused by the sudden change of groundwater level during heavy rain period can be accurately captured.

[0008] Dynamic subdomain discretization technology is used in sandy soil areas. Its elastic modulus correction model can update the attenuation effect of changes in porosity and water content on soil stiffness in real time. In clay areas, a hybrid registration algorithm is used to separate rigid displacement from plastic creep, solving the coupling problem of consolidation settlement and structural deformation that is difficult to distinguish using traditional methods, and keeping the cumulative error of long-term monitoring within 0.3 mm. Through the time-frequency domain feature extraction technology, the constructed spatiotemporal correlation matrix can reveal the feedback relationship between the displacement field and geological parameters (such as cohesion and internal friction angle). This method can warn of potential instability risks. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 It is a flow chart of an intelligent monitoring method for displacement of a stormwater resource interception and storage structure provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0010] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying 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 to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0011] like Figure 1 As shown, an embodiment of the present invention provides a method for intelligently monitoring the displacement of a rainwater resource interception and storage structure, the method comprising the following steps: Step 1: Based on the initial displacement field constraints obtained by point cloud registration, a three-dimensional spatial distribution model is constructed by coupling multi-source geological data through particle swarm optimization; the sandy soil area is discretized into dynamic subdomain grids, and parameter values ​​are assigned to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid; based on the parameterized grid, an elastic modulus correction model is established; based on the elastic modulus correction model, the deformation field of the particle swarm is iterated to obtain the adjusted deformation field; and the displacement correction factor is calculated based on the adjusted deformation field; Step 2, obtaining the second source point cloud and the second target point cloud of the clay area, and establishing a spatial distribution map; 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 map and the global rigid transformation matrix to obtain the registered point cloud; calculating the corrected displacement field sequence based on the registered point cloud; Step 3: Perform time-frequency domain analysis on the displacement correction factor and the corrected displacement field sequence, extract characteristic frequencies and damping ratios, construct a spatiotemporal correlation matrix, and generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters.

[0012] In the embodiment of the present invention, step 1, the sandy soil is divided into 5-8 layers of variable grids, each layer is assigned dynamic porosity and moisture content parameters, and the real-time solution of the permeability-effective stress coupling response during the rainstorm period is realized, and the simulation accuracy is improved by 45%. The iterative particle swarm algorithm is introduced to optimize the soil stiffness attenuation coefficient, and the deformation underestimation problem caused by the constant modulus value in the traditional method is solved, so that the correlation coefficient between the predicted settlement 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 the explanation of deformation mechanism. Compared with the empirical formula method, the creep prediction accuracy is improved by 73%; the rigid matrix corrects the overall translation / rotation error, and the plastic gradient function compensates for the local shear band deformation. Double registration improves the smoothness of the displacement field in the clay region by three orders of magnitude; through the space-time filtering technology, environmental noise such as temperature gradient and vegetation growth is eliminated to obtain a pure deformation signal with a signal-to-noise ratio higher than 40dB, which meets the needs of millimeter-level deformation monitoring; the inherent period of the displacement field (such as daily period and rainy season period) is identified, and the load response spectrum is established, which can 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 is greater than 0.05, an early warning is triggered, which reduces the false alarm rate by 82% compared with the traditional threshold method, and reveals the space-time coupling law of displacement hotspots and geological parameters (such as cohesion and compression modulus), providing a decision-making basis for targeted reinforcement and optimizing maintenance costs by about 40%.

[0013] In a preferred embodiment of the present invention, 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, may include: 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 3D coordinates of the measuring points, including: using RieglVZ-400 3D 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 a point density of ≥300pts / m 2 , vertical accuracy ±3mm, horizontal accuracy ±5mm, deploy Geoprobe7822DT geological radar system, arrange 5 measuring lines along the main sliding direction of the slope, with a measuring point spacing of 2m, and use Guelph permeameter to conduct in-situ permeability test to obtain the permeability coefficient K (accuracy ±0.5×10 -5 m / s), use a nuclear density meter to measure the porosity n (accuracy ±1.5%), and record the three-dimensional coordinates (X, Y, Z) of each measuring point through RTK-GPS measurement, with a plane accuracy of ±1cm and an elevation accuracy of ±2cm.

[0014] Step 12: align the denoised source point cloud and target point cloud and calculate the initial displacement field, including: applying statistical outlier removal filtering, setting the neighborhood point threshold to 50, the standard deviation multiple to 1.5, using the moving least squares (MLS) method for surface smoothing, the smoothing window size to 0.3m, using the improved ICP algorithm (Generalized-ICP), combined with the normal vector matching strategy, setting the maximum number of iterations to 50, and the convergence threshold to 10. -5 m, use KD-Tree to accelerate the nearest neighbor search and calculate the initial displacement field: △X = Xt-Xs, △Y = Yt-Ys, △Z = Zt-Zs, where (Xs, Ys, Zs) are the coordinates of the point in the source point cloud, and (Xt, Yt, Zt) are the coordinates of the corresponding point in the target point cloud.

[0015] Step 13: convert the geological exploration data of each measuring point into a particle, and the particle attributes include permeability, porosity and spatial coordinates; generate a particle group in three-dimensional space, and simulate the spatial correlation of geological parameters through virtual "force", which includes: The geological exploration data of each measuring point is converted into a particle. The particle attributes include permeability, porosity and spatial coordinates. Each particle represents a geological exploration point. The attribute vector ; Initialize the particle position to the measured coordinates, add ±5% random perturbations to the permeability coefficient and porosity, which can increase the diversity of the particle group and avoid falling into the local optimal solution; Virtual force modeling: Gravitational term: , simulating the spatial continuity of geological parameters, where is the gravitational coefficient, and is the “mass” of the particle (which can be defined according to actual conditions), is the distance between particles, is the unit direction vector between particles; Repulsion term: , to prevent excessive aggregation of particles, where k is the repulsion coefficient, is the equilibrium distance between particles; Viscosity term: , simulating the damping effect of geological processes, is the viscosity coefficient, is the velocity of the particle; Parameter settings, number of particles N = 1000 (dynamically adjusted according to the density of measurement points), virtual force coefficient: G = 0.01, k = 100, β = 0.8, neighbor particle search radius .

[0016] Step 14, calculate the evaluation value of the geological parameters of the current particle position, update the particle speed and position, and stop the optimization when the particle group converges to obtain the converged particle group, which specifically includes: updating the particle speed and position according to the current position, speed and virtual force of the particle. The specific update formula is as follows: Speed ​​Update: ; Location Updates: ; in, is the inertia weight, and is the learning factor, and is a random number, is the best historical position of the particle, is the global optimal position, is the sum of the virtual forces acting on the particles; when the particle swarm meets the convergence condition, the optimization is stopped. The convergence condition can be that the change in the particle position is less than a certain threshold, and the converged particle swarm is obtained.

[0017] Step 15, based on the converged particle swarm, generates a continuous three-dimensional spatial distribution model of permeability coefficient and porosity, specifically including: Obtain the converged particle swarm data 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 interpolation algorithm input, such as taking 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 value of the unknown point through the data of the known point, and use the semivariogram function to describe the spatial correlation between the data. The goal of Kriging interpolation is to find a set of weight coefficients that minimize the estimation error of the interpolation result.

[0018] Calculate the spatial distance between each pair of particles in the particle swarm , for each pair of particles, calculate half the difference in their properties (permeability or porosity), i.e. ,in and Particles and particles Attribute value of The calculated distance and semivariance , select a suitable semivariogram model (such as a 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, the Kriging equations are constructed based on the semivariogram model and the positions of the known points. The general form of the Kriging equations is: ; in, It is a known point and The semivariance between is the point to be interpolated With known points The semivariance between 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 and obtain the weight coefficients and Lagrange multipliers .

[0019] For each interpolation point, according to the weight coefficient obtained and the property values ​​(permeability or porosity) of the known points to calculate the interpolation results , in three-dimensional space, the area to be interpolated is traversed according to a certain grid spacing (such as set according to the required model resolution), and the above interpolation calculation is performed on each grid point to obtain the permeability coefficient and porosity value of the point.

[0020] The permeability coefficient and porosity values ​​in the interpolated three-dimensional space are organized into a data format suitable for visualization and subsequent analysis, such as voxel data (Voxeldata) or grid data (Griddata). Professional geographic information system (GIS) software, geological modeling software or programming language (such as Python's Mayavi, VTK and other libraries) are used to generate a continuous three-dimensional spatial distribution model from the organized data. The model can display the distribution of permeability coefficient and porosity in space in the form of three-dimensional graphics.

[0021] In the embodiment of the present invention, measuring points are deployed in a 20m×20m grid to obtain parameters such as permeability coefficient K and porosity n, and a three-dimensional geological database containing 15-20 layers of geological interfaces is constructed, which is 8 times denser than traditional drilling data. A millimeter-level precision three-dimensional laser scanning system is used to obtain the source point cloud (initial state) and the target point cloud (deformed state), with a point density of 500pts / m 2 , completely retaining the micro-topography features; Step 12, use the cloth simulation filtering algorithm to eliminate vegetation jitter noise, improve the signal-to-noise ratio by 12dB, and achieve global optimal matching of source-target point clouds through iterative nearest point search. The initial displacement field accuracy reaches sub-millimeter level, which is 60% more efficient than the traditional optical flow method. Step 13, convert each measurement point data into particles, assign physical properties such as mass and charge, and the virtual interaction force between particles matches the measured spatial correlation function by 92%.3 Generate 5000-8000 particle groups in space, simulate the spatial continuity of geological parameters through repulsion / attraction balance, and avoid the over-smoothing phenomenon of traditional interpolation methods. Step 14, integrate the variation function in geostatistics and the point cloud displacement constraint, construct a multi-objective optimization function, so that the permeability field and the measured groundwater flow field are consistent with 95%, adopt a linear decreasing inertia weight strategy, focus on global search (w=0.9) in the early stage, and strengthen local development (w=0.4) in the later stage, and the convergence speed is 40% higher than that of standard PSO. Step 15, select the Multi-Quadric function to reconstruct the surface of the converged particle group, generate a continuous parameter model with a resolution of 0.5m, which is 7 times more efficient than the Kriging interpolation method, and simultaneously generate the three-dimensional distribution of permeability, porosity and effective stress field, so as to realize the data basis of seepage-deformation coupling analysis.

[0022] In a preferred embodiment of the present invention, calculating the estimated value of the geological parameter of the current particle position includes: Based on the porosity of the current particles and the effective particle size of the soil particles, the theoretical permeability coefficient of the particle position is calculated, specifically including: The theoretical permeability coefficient K is calculated based on the porosity n and the effective particle size d10 of the soil particles. The theoretical empirical formula is the Kozeny-Carman equation: ; in, is the effective particle size of soil particles (unit: m), usually the particle size corresponding to a cumulative mass percentage of 10%, and n is the porosity (dimensionless), ranging from 0 to 1.

[0023] Compare the measured permeability coefficient of the current particle with the theoretical permeability coefficient, calculate the absolute value of the relative deviation, and obtain the physical constraint value, specifically including: obtaining the measured permeability coefficient of the current particle , the theoretical permeability coefficient is obtained from the above steps ,Will and Substitution , calculate the absolute value of the relative deviation , this value is the physical constraint value.

[0024] Taking the current particle as the center, according to the preset neighborhood search radius, determine the set of neighboring particles around it, and calculate the standard deviation of the permeability coefficient and porosity within the neighborhood search radius; according to the standard deviation of the permeability coefficient and porosity, calculate the spatial constraint value, specifically including: assuming the position of the current particle to be (x0, y0, z0), the preset neighborhood search radius to be r, traversing all particles, for particle i, its position is (xi, yi, zi), calculate the Euclidean distance di between particle i and the current particle, if di≤r, particle i belongs to the neighborhood particle set of the current particle; assuming the value of the permeability coefficient in the neighborhood particle set to be K1, K2, ..., Km, calculate the average value of the permeability coefficient, and calculate the standard deviation of the permeability coefficient based on the average value; assuming the value of the porosity in the neighborhood particle set to be n1, n2, ..., nm, 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: ; in, and They are the weights of the standard deviation of permeability coefficient and porosity, respectively, and can be set according to actual conditions, for example , .

[0025] The physical constraint value and the spatial constraint value are weighted and summed to obtain the evaluation value, which specifically includes: The calculation formula of the evaluation value E is: ; in: is the weight of the physical constraint value (dimensionless), ranging from 0 to 1; is the weight of the spatial constraint value (dimensionless), and the sum of the two weights is 1.

[0026] In the embodiment of the present invention, the method not only takes into account the geological parameters of the current particle itself, but also takes into account the spatial variability through neighborhood analysis, so that the evaluation result is more comprehensive and reliable. By comparing the measured permeability coefficient with the theoretical permeability coefficient, it is ensured that the evaluation result conforms to the physical laws, thereby improving the scientificity and accuracy of the evaluation. Parameters such as weights and neighborhood search radius in the method can be adjusted according to actual conditions, so that the evaluation method has greater flexibility and applicability, and can be applied to different geological conditions and evaluation needs. The evaluation value provides a quantitative basis for geological engineering design and decision-making, which helps to optimize the design scheme and improve engineering efficiency.

[0027] In a preferred embodiment of the present invention, the sandy soil area is discretized into dynamic subdomain grids, and parameter values ​​are assigned to each subdomain based on a three-dimensional spatial distribution model to form a parameterized grid, including: The three-dimensional spatial distribution model is converted into three-dimensional grid data, which includes: determining the size of the grid (such as the side length is dx, dy, dz) and the number of rows, columns, and layers of the grid according to the scope and accuracy requirements of the study area. For example, if the study area ranges in the x direction , the y direction range is , the z-direction range is , the grid side lengths are dx, dy, dz, then the number of grids in the x direction is , the number of grids in the y direction is , the number of grids in the z direction is ; For each point in the three-dimensional spatial distribution model , calculate the grid cell index where it is located , the calculation method is: ; ; ; ; in, Indicates rounding down. The geological parameter values ​​(such as permeability, porosity, etc.) are assigned to the corresponding grid cells. , if there are multiple points in a grid cell, the parameter value of the grid cell can be determined by using methods such as average value, maximum value, and minimum value.

[0028] According to the filtered 3D grid data and soil interface geometry information, the subdomains are divided step by step on the 2D plane to obtain dynamic subdomain grids, including: Filter the three-dimensional grid data to remove noise and outliers. Mean filtering, median filtering and other methods can be used. For example, for each grid cell, the average or median value of the grid cell parameter values ​​in a certain neighborhood around it is taken as the filtered value of the grid cell. Extract the geometric information of the soil interface from the three-dimensional spatial distribution model, such as the equation, vertices and boundaries of the soil interface. The edge detection algorithm in image processing (such as the Canny algorithm) or the interface recognition algorithm in geological modeling can be used to extract the soil interface. Select a suitable two-dimensional plane (such as a horizontal plane or a vertical section) and divide the subdomains on the plane according to the soil interface information. The data structure such as quadtree or octree can be used for step-by-step division. For example, for quadtree partitioning, the two-dimensional plane is first divided into four sub-areas, and then each sub-area is judged. If the geological conditions in the sub-areas are quite different, the sub-area is further divided into four smaller sub-areas until the division accuracy requirement is met. According to the sub-domain division results on the two-dimensional plane, combined with the three-dimensional grid data, a dynamic sub-domain grid is generated. Each sub-domain corresponds to a part of the grid unit in the three-dimensional space, and the boundary of the sub-domain can be determined based on the division boundary on the two-dimensional plane and the vertical direction of the three-dimensional grid data.

[0029] According to the permeability and porosity values ​​of the dynamic subdomain grid and the corresponding grid cells, each subdomain is analyzed to calculate the average and maximum values ​​of the permeability, the average and minimum values ​​of the porosity, and the coefficient of variation of the permeability and porosity within the subdomain, including: For each subdomain, traverse the corresponding grid cells, extract the permeability and porosity values, and calculate the average permeability value according to the number of grid cells in the subdomain and the permeability values ​​of the grid cells; calculate the average porosity according to the porosity values ​​of the grid cells, calculate the maximum permeability value and the minimum porosity; calculate the ratio of the standard deviation of the permeability coefficient to the average permeability coefficient to obtain the permeability variation coefficient; calculate the ratio of the standard deviation of the porosity to the average porosity to obtain the porosity variation coefficient.

[0030] The average value, maximum value, minimum value and coefficient of variation of the subdomain are combined into a parameter vector as the geological attribute of the subdomain, including: for each subdomain, the average value of its permeability , maximum value The average porosity , minimum , coefficient of variation of permeability and the coefficient of variation of porosity Combined into a parameter vector , store the parameter vector for each subdomain in a data structure such as a list or dictionary.

[0031] The parameter vector is normalized by minimum-maximum normalization to obtain a normalized parameter vector.

[0032] The normalized parameter vector is used as a 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, the subdomain to which each grid cell belongs is determined according to the division result of the dynamic subdomain grid, and the normalized parameter vector of the corresponding subdomain is used as a geological attribute parameter and attached to each grid cell in the subdomain. A data structure (such as a list or array) is created 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 a dynamic subdomain grid, and parameter values ​​are assigned to each subdomain to form a parameterized grid.

[0033] In an 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; sub-domain division combined with the geometric information of the soil layer interface can better reflect the geological structure characteristics of the sandy soil area. The physical and mechanical properties of different soil layers may be different. Dividing the sub-domains according to the soil layer interface can make the geological conditions in each sub-domain relatively uniform, thereby improving the accuracy of subsequent parameter analysis. The dynamic subdomain grid is obtained by a step-by-step division method. The size and shape of the subdomain can be flexibly adjusted 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 characteristics; 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 the parameters in the subdomain, and the maximum and minimum values ​​reveal the extreme conditions of the parameters, which is helpful to evaluate the geological stability and carrying capacity of the subdomain. The calculated coefficient of variation can measure the degree of discreteness of the permeability and porosity in 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 in the subdomain and the more complex the geological conditions; the smaller the coefficient of variation, the relatively uniform parameters in the subdomain and the relatively stable geological conditions, which is helpful to identify areas with abnormal geological conditions. Different statistical indicators (such as mean, maximum, minimum, and coefficient of variation) have different dimensions and value ranges. Normalization can eliminate the impact 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 geological attribute parameters with the vertex coordinates of the grid unit, realizing the precise association between geological information and spatial position. The vertex coordinates of the grid unit can accurately locate the position of the geological attribute parameters in three-dimensional space.

[0034] In a preferred embodiment of the present invention, an elastic modulus correction model is established according to the parameterized grid; according to the elastic modulus correction model, the deformation field of the particle group is iterated to obtain an adjusted deformation field, including: Determine the joint distribution constraint relationship between porosity and permeability, specifically including: collecting experimental data of permeability under different porosities (such as sandstone, carbonate rock, etc.), using regression analysis (such as linear regression, polynomial regression) to fit the relationship curve between porosity n and permeability k. For example, for sandstone, the relationship can be obtained: ,in, and For the fitting constant, the theoretical model is verified: The rationality of experimental data is verified by combining theoretical formulas (such as the Kozeny-Carman equation). The Kozeny-Carman equation is expressed as: ;in, is the particle diameter, For tortuosity, by comparing experimental data with theoretical predictions, the model parameters were adjusted to improve accuracy.

[0035] Obtain the measured data of elastic modulus and porosity of different lithology samples, analyze their corresponding relationship, and establish the lithology function relationship between elastic modulus and porosity for different lithology, 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 functional relationship between elastic modulus and porosity. For example, for sandstone, the relationship can be obtained: ;in, is the elastic modulus when the porosity is 0, is a constant.

[0036] According to the joint distribution constraint relationship and the lithology function relationship, the elastic modulus correction model of the joint constraint of permeability coefficient and porosity is calculated, which specifically includes: combining the porosity-permeability coefficient relationship (such as ) and elastic modulus-porosity relationship (e.g. ), an elastic modulus correction model is established, wherein the calculation formula of the elastic modulus correction model is: ,in, represents the modified elastic modulus; is the elastic modulus when the porosity is 0; is the constant of elastic modulus changing with porosity; is the porosity; is the constant in the relationship between porosity and permeability; is the base of natural logarithms; Represents the linear adjustment parameter in the deformation field.

[0037] According to the elastic modulus correction model, the permeability coefficient and porosity data of the parameterized grid, the static elastic modulus value of each grid unit is calculated to form the initial elastic modulus field, which includes: according to the porosity n and permeability coefficient k data of the parameterized grid, the correction model is applied to calculate the elastic modulus of each unit .

[0038] The elastic modulus values ​​of all grid cells are combined into a three-dimensional spatial distribution model, and the simulated deformation field is calculated, including: discretizing the geological body using octree (Octree) or tetrahedral grid (TEN), defining the grid cell size, mapping the elastic modulus value of each grid cell (calculated by the modified model) to the corresponding grid node, and using Kriging interpolation or inverse distance weighted interpolation to generate a continuous field for irregularly distributed porosity / permeability data. Formula example (inverse distance weighted interpolation): ;in, is the elastic modulus of the point to be interpolated, is the elastic modulus of the known point, is the power parameter (usually 2), N represents the number of known data points involved in the interpolation calculation, Represents the distance between the point to be interpolated and the i-th known point.

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

[0040] Read data from files storing three-dimensional elastic modulus field data (such as CSV files, which contain the elastic modulus values ​​corresponding to each grid node or unit) and map them to the corresponding grid nodes or units. For irregularly distributed data, use Kriging interpolation or inverse distance weighted interpolation 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.

[0041] 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 unit, and update the material properties. Ensure that the elastic modulus value of each node or unit is consistent with the actual geological conditions; determine the nodes or surfaces that need to be subject to displacement constraints based on the input boundary condition information, and apply the corresponding displacement constraints (such as fixing the displacement in the X, Y, and Z directions to 0) to the specified nodes or surfaces.

[0042] Determine the nodes or surfaces where tectonic stress needs to be applied, and apply tectonic stress to the corresponding nodes or surfaces according to the input stress value and direction. Enable the gravity load function and set the direction and magnitude of gravity acceleration (usually -Z direction, acceleration is 9.81m / s 2 ), ANSYS automatically calculates gravity loads based on the density of the material.

[0043] According to the scale and complexity of the model, select a suitable solver, such as a sparse solver; set the convergence accuracy of the solution, the maximum number of iterations and other parameters, for example, set the force convergence tolerance to 1%-0.1%, the maximum number of iterations to 50-100 times, and submit the configured model and boundary conditions to the ANSYS solver for calculation. During the solution process, monitor the number of iterations, residual ratio and other information in real time. If the residual ratio continues to fail to converge or exceeds the preset maximum number of iterations, stop the calculation and output an error message. If an abnormal situation occurs during the solution process, such as mesh distortion, missing material properties, etc., automatically diagnose the error and try to take appropriate measures to repair it, such as adjusting the mesh quality, checking material properties, etc.

[0044] After the solution is completed, the displacement data of each grid node is extracted from the calculation results of ANSYS, including the displacement components in the X, Y, and Z directions. The extracted displacement data is post-processed, such as calculating the total displacement, and the processed data is stored in the specified file.

[0045] When applied specifically, the simulation process may include: defining material parameters: Poisson's ratio ,density Read from input file, elastic modulus ; Displacement boundary: Fixed base node displacement (such as ), the constraint equation is: ; Stress boundary: Apply structural stress at the edge of the model , such as horizontal stress . Gravity Load: Enable gravity acceleration , body force Solver configuration and numerical calculation: Select a sparse solver according to the model size and set the convergence conditions: Force convergence tolerance: (relative residual), Represents the residual vector after the kth iteration, indicating the current approximate solution U k The unbalanced amount after substituting into the equation; Represents the load vector The L2 norm of the node; the maximum number of iterations: 50~100. 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 node displacement vector (the deformation field to be solved). Read the three-dimensional displacement components of each node from the solution results , forming the displacement vector: ; Where N represents the total number of grid nodes, and the total displacement is calculated for each node , forming displacement cloud map data, storing the displacement data as CSV / EXCEL files according to node numbers, coordinates, and component values, and finally obtaining the simulated deformation field.

[0046] Compare the simulated deformation field with the actual observation 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 pore pressure field changes, and update the elastic modulus field through the pore pressure changes. When the deformation field adjustment amount 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, including: Get calculation results from ANSYS, including node coordinates , 3D simulation displacement (The result of the kth iteration is recorded as ); displacement data obtained through point cloud registration (sandy soil) or mixed registration (clay), including measurement point coordinates and displacement .

[0047] Convert the simulation and observation data into the same coordinate system (such as the local coordinate system of the project), unify the units (such as meters), and ensure the consistency of the coordinate origin and scale; calculate the coordinate deviation of the same-name points, requiring the maximum deviation to be less than 1 / 10 of the grid unit size (for example, when the grid accuracy is 0.5m, the deviation must be less than 0.05m).

[0048] For unstructured observation point clouds (such as clay area point clouds), inverse distance weighted interpolation or Kriging interpolation is used to map the observed displacements to the simulation grid nodes to generate Isomorphic observed displacement fields .

[0049] For each mesh node i, calculate the displacement error in three directions: ; ; ; Where k is the iteration number, i is the node number, and reflects the displacement deviation between simulation and observation in the x, y or z direction (unit: m). Total displacement absolute error: ; Represents the spatial sum error between the single-point simulation and the observed displacement, which is used to intuitively judge the local matching accuracy (unit: m).

[0050] Subdomain level error statistics (sandy soil area): Group calculations by subdomain of the parameterized grid (Total M subdomains), count the number of The error of each node includes: Mean absolute error: (Reflects the overall error level of the subdomain).

[0051] Error standard deviation: (Measures the degree of error dispersion within a subdomain).

[0052] Coefficient of Variation: (dimensionless, indicates significant error fluctuation).

[0053] Geological parameter correlation: Extract high error subdomains ( , is the global error standard deviation) , porosity , analyze whether the error is related to parameter anomalies (such as The simulated value is 30% lower than the measured value).

[0054] Point cloud cluster difference analysis (clay area): Divide the point cloud into high cohesion areas according to the cohesion map and low cohesion zone , comparing the mean values ​​of the combined displacement errors of the two types of regions, verifying the effectiveness of the plastic deformation constraint (theoretically, the error in the high cohesion area should be smaller), where: represents the average cohesion, Indicates cohesion; Standard deviation of cohesion.

[0055] Annotation on 3D geological model , mark the "abnormal area" where the error exceeds the threshold (such as 2 times the global standard deviation), superimpose the permeability coefficient contour surface or cohesion contour line, and check whether the error anomaly coincides with the geological interface (such as weak interlayer); draw the combined displacement curve of simulation and observation along the key profile of the structure (such as the central axis of the dam foundation), mark the position where the difference exceeds the allowable error of the project (such as 10mm), and locate the specific deformation out of control area.

[0056] With the goal of minimizing the global sum of squared errors, define: ; in, is the total number of grid nodes, that is, the number of nodes after discretization in the calculation area; The smaller the value, the better the match between simulation and 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 modulus or the 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.

[0057] The sensitivity of parameters to errors is calculated by the finite difference method, such as: ; For example, fixed and , perturbation get The rate of change of ,in, represents the first Row, No. The elements of the column represent The total displacement error of the nodes is Parameters The partial derivative of Represents the model parameters to be optimized.

[0058] Parameter increment calculation: ,in, Represents the total displacement error vector, with dimension N×1; The Gram matrix (symmetric positive definite matrix) representing the sensitivity matrix, is the pseudo-inverse matrix used to solve the linear least squares problem.

[0059] Solve the parameter increment through matrix operation, so that The fastest decrease; for each node, the relaxation factor λ (0.1~0.5, to avoid over-adjustment) is used to adjust the displacement: (similar to 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 kth iteration; Indicates In the iteration, the updated x-direction displacement (the same applies to the y and z directions)

[0060] For subdomains where the mean error exceeds the standard , a rigid displacement correction is applied along the main deformation direction n: ,in, represents the mth subdomain; Indicates subdomain The average displacement error modulus of all nodes in the subdomain (scalar, reflecting the size of the error); n represents the average direction unit vector of the displacement vector in the subdomain (vector, pointing to the main deformation direction, obtained by averaging the normalized displacement vectors of each node in the subdomain); Represents the rigid displacement correction amount (vector) of the entire subdomain, applied along the main deformation direction n, and used to correct the systematic deviation of the local area.

[0061] Pore ​​pressure field inversion and elastic modulus update, pore pressure field inversion (based on Darcy's law); Seepage control equation: ; in, represents the permeability coefficient (m²), from the parameterized grid or the modified model; Indicates the dynamic viscosity of the fluid (Pa・s, such as water, 10 −3Pa·s); p represents the pore pressure Pa, reflecting the pressure of pore water on the soil skeleton; represents the water storage rate (1 / m), which is related to the compressibility of soil and fluid.

[0062] The porosity is dynamically updated by the adjusted deformation field Calculate the volumetric strain (i.e., change in porosity): ; in, Indicates soil expansion (increase in porosity), indicates compression (porosity reduction), represents the porosity change; Indicates In the iteration The three-dimensional displacement field of node i in a subdomain contains three components: , and Respectively represent the nodes in , Displacement in the and directions (unit: m); represents the divergence of the displacement field, i.e., the volume strain. The permeability coefficient is updated: Based on the Kozeny-Carman equation, the permeability coefficient is positively correlated with the cubic porosity: ; in, Indicates The permeability coefficient of the iteration (unit: m 2 ); Indicates The permeability coefficient after the iteration update (unit: m 2 ); Indicates The porosity of the iteration; Indicates The porosity of the iteration reflects the porosity change after soil deformation; for example, a 10% increase in porosity will increase the permeability by about 33% ( ).

[0063] For each grid cell, substitute the updated porosity and model parameters ; ; Among them, Em(k+1) represents the elastic modulus (unit: Pa) after the k+1th iteration update, reflecting the stiffness of the soil; E0 represents the initial elastic modulus (unit: Pa), that is, the theoretical elastic modulus when the porosity is zero; express The increment of and Respectively and The increment of represents the updated porosity; Example: If , the effect of porosity on elastic modulus is aggravated, and the high porosity area Further reduce the updated Assigned to the corresponding grid nodes / elements to form a new elastic modulus field, which serves as the material property input for the next round of ANSYS calculations.

[0064] Calculate the global displacement adjustment rate to determine whether the deformation field is stable: ; is a preset threshold (such as 0.5%), indicating that the adjustment amount is considered small enough when the overall change of the deformation field is less than 0.5%. It represents the global displacement adjustment rate, which measures the relative change of the deformation field between two adjacent iterations; Indicates Simulated displacement of iteration (unit: m); difference threshold of consecutive iterations: Compare the root mean square error (RMSE) of two consecutive iterations to ensure that the error fluctuation is within a controllable range: ; in, represents the root mean square error (unit: m), which measures the overall deviation of the simulated displacement from the observed displacement; Indicates the RMSE fluctuation threshold (for example, the preset value is 0.1mm); if the number of iterations exceeds the preset maximum value (for example, 20 times), it is forced to terminate even if it has not converged to avoid wasting computing resources. Node displacement field that meets the convergence conditions , contains the three-dimensional displacement of each node , used for subsequent displacement correction and structural safety assessment.

[0065] Updated elastic modulus field: 3D distribution model including the latest geological parameters , record the RMSE and parameter adjustment of each round , changes in pore pressure peak values, etc., are used to trace the model optimization process and provide a reference for parameter adjustment for subsequent similar projects.

[0066] In the embodiment of the present invention, porosity and permeability are important parameters for describing the characteristics of geological materials. There is an intrinsic connection between them. Determining the joint distribution constraint relationship can accurately reflect the coupling relationship between the internal pore structure of the geological body and the fluid flow characteristics, which is helpful to simulate the mechanical and hydraulic behavior of the geological body more realistically. In the subsequent elastic modulus correction model and deformation field calculation, considering the joint distribution constraint relationship between porosity and permeability can make the model more consistent with the actual geological conditions and improve the accuracy and reliability of the simulation results. Geological materials of different lithologies have different physical and mechanical properties, and there are also differences in the relationship between their elastic modulus and porosity. Establishing a lithology function relationship can fully consider this difference and make the model more targeted and applicable. The measured data provides a reliable basis for the establishment of the lithology function relationship, ensures the accuracy and reliability of the function relationship, and through the analysis of a large amount of measured data, the inherent law between the elastic modulus and the porosity can be revealed. Combining the joint distribution constraint relationship between porosity and permeability and the lithology function relationship between elastic modulus and porosity, the hydraulic and mechanical properties of the geological body are comprehensively considered, making the elastic modulus correction model more comprehensive and accurate. The revised model can automatically adjust the calculation of the elastic modulus according to different geological conditions (such as porosity, permeability and lithology), improve the adaptability of the model to different geological environments, and make the simulation results closer to the actual situation. By revising the model, the errors caused by not considering the coupling relationship of multiple factors in the traditional model can be eliminated, the accuracy and reliability of the model can be improved, and a more accurate basis can be provided for the subsequent deformation field calculation. The parameterized grid discretizes the study area into multiple grid units, and the elastic modulus value of each grid unit is calculated to achieve the discretization distribution of the elastic modulus in space.

[0067] In a preferred embodiment of the present invention, the displacement correction factor is calculated based on the adjusted deformation field, including: According to the geological stratification information of the sandy soil area, the deformation field data is divided into several horizontal layers, including: extracting the stratification information of the sandy soil area from the geological survey report, the borehole column chart or the three-dimensional geological model, including the top elevation, bottom elevation, lithology description (such as silt, fine sand, medium sand), sedimentary age, etc. of each soil layer; dividing by elevation interval (such as every 5 meters as a horizontal layer) or by geological unit (such as natural soil layers divided according to sedimentary cycles) to ensure that the geological properties (such as porosity and permeability coefficient) in each layer are similar; extracting the elevation z of each point for the adjusted deformation field data, including the three-dimensional coordinates (x, y, z) and displacement vector (ux, uy, uz) of each grid node; assigning the points in the deformation field to the corresponding horizontal layers one by one according to the elevation range of the geological stratification. For example, if the elevation range of a layer is 10m≤z<20m, all points that meet this condition are classified as this layer; For each layer, the displacement vectors of all points in the layer are calculated, and the characteristic displacement of the layer is calculated. The characteristic displacement includes the maximum displacement, the minimum displacement and the average displacement. Specifically, for all points in each layer, the displacement components are extracted, the modulus (total displacement) of the resultant displacement vector is calculated, the total displacement of all points is traversed layer by layer, and the maximum and minimum values ​​are taken; the arithmetic mean of the total displacement of all points in the layer is calculated.

[0068] Obtain the sand state parameters of the sandy soil area and spatially match the sand state parameters with the deformation field data so that the displacement vector of each point has its corresponding sand state parameters, including: Obtain sand state parameters from geotechnical test data (such as standard penetration test, moisture content test, particle analysis), including: density (relative density, dimensionless), moisture 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: Using geological drilling data, the sand state parameters of adjacent boreholes are obtained. Spatial interpolation methods (such as inverse distance weighted interpolation and Kriging interpolation) are used to calculate the sand state parameters of the point according to the parameter values ​​of adjacent boreholes. The displacement vector of each point and the corresponding sand state parameters are stored in structured data (such as a table or database) to ensure a one-to-one correspondence.

[0069] The displacement correction coefficient of each layer is calculated using the characteristic displacement and sand state parameters, including: , average water content ) and characteristic displacement , construct the correction coefficient calculation model, including: According to the mechanical properties of sand, the negative correlation between the correction coefficient and density is set, for example: ; Where k is the correction coefficient sensitivity parameter (calibrated by test or historical data), is the average density within the layer (value range 0≤ ); Input the characteristic displacement of each layer (such as the average displacement ) and sand state parameters (such as average density , average water content ), substitute it into the correction model to calculate the displacement correction coefficient of this layer ; Example: If the average density of a layer of sand is =0.6 (medium density state), the correction coefficient is calculated according to the empirical formula =0.9, indicating that the simulated displacement of this layer needs to be multiplied by 0.9 to correct it to be closer to the actual displacement.

[0070] According to the adjusted deformation field and displacement correction coefficient, the displacement correction factor is calculated, which specifically includes: The correction factor can be expressed as: ,in, 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 density anomaly points); for all points in each layer, the correction coefficient of the layer is As a basic correction factor. If point-level refinement is required, it can be based on the single-point sand state parameters (such as the density of a certain point). 10% lower than the layer average). To make fine adjustments: ,in is the point-level correction weight (determined by sensitivity analysis). Multiplying with the original displacement vector gives the corrected displacement: , and finally a three-dimensional deformation field including the correction factor is formed.

[0071] In the embodiment of the present invention, the sandy soil area is divided into several horizontal layers, and the displacement characteristics of different soil layers can be considered more finely. Since the geological conditions of each layer are relatively uniform, the displacement vectors of all points in the layer can be calculated more accurately after stratification, thereby improving the calculation accuracy of the overall displacement field. By calculating the maximum displacement, the minimum displacement and the average displacement, the displacement distribution of each layer can be fully grasped. These characteristic displacements not only reflect the extreme value and average level of displacement, but also help identify potential risk areas or abnormal displacements. The sand state parameters (such as density, porosity, water content, etc.) are spatially matched with the deformation field data to ensure that the displacement vector of each point has its corresponding sand state parameters. This matching makes the displacement calculation more in line with the actual situation, because the physical and mechanical properties of 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. Using the characteristic displacement and the sand state parameters to calculate the displacement correction coefficient, the influence of different factors on the displacement can be quantified. This comprehensive consideration makes the corrected deformation field more in line with the actual situation and enhances the reliability of the model.

[0072] In a preferred embodiment of the present invention, a second source point cloud and a second target point cloud of the clay region are obtained, and a spatial distribution map is established; based on the second source point cloud and the second target point cloud, a global rigid transformation matrix is ​​calculated, including: Collect original point cloud data in the clay area, denoise the original point cloud to obtain denoised point cloud data; segment the second source point cloud and the second target point cloud corresponding to the clay area through the regional growing algorithm according to the denoised point cloud data, and downsample the segmented point cloud to obtain the second source point cloud and the second target point cloud, which specifically includes: using a three-dimensional laser scanner to collect data according to a certain scanning strategy to ensure that the entire clay area is covered, and pay attention to setting appropriate scanning parameters (such as resolution, scanning angle, etc.) during the collection process; calculate the average distance from each point to its neighboring points, set a threshold value according to the statistical distribution of the average distance, determine the points with too large distance as noise points and remove them, set a radius range with each point as the center, count the number of neighboring points within the radius, and if the number of neighboring points is less than a certain threshold, determine the point as a noise point and remove it; Some representative points are selected manually or automatically as seed points according to the features of the point cloud (such as curvature, normal direction, etc.). For each seed point, its neighborhood points are searched, and the preset growth criteria (such as normal angle, distance between points, etc.) are used to determine whether the neighborhood points belong to the same area. If the criteria are met, the neighborhood point is added to the current area, and the growth is continued with the neighborhood point as the center until there are no neighborhood points that meet the conditions. Through multiple growths, the point cloud is divided into different areas, from which the second source point cloud and the second target point cloud corresponding to the clay area are screened out, and the point cloud space is divided into several small voxel grids. For each point in the voxel grid, its center of mass is calculated, and all points in the voxel grid are replaced by the center of mass point, thereby reducing the number of point clouds.

[0073] Obtain clay samples at different depths and the corresponding cohesion parameters; establish a three-dimensional distribution model of porosity and water content based on the clay samples and the corresponding cohesion parameters; establish a spatial distribution map of cohesion using the collaborative Kriging spatial interpolation method based on the three-dimensional distribution model of porosity and water content, specifically including: collect undisturbed soil samples by drilling at different positions and depths (e.g., every 2 meters) in the clay area, each sample weighing about 200-500 grams, immediately seal them in a moisture-proof bag, and record the sample collection coordinates (X, Y, Z axis position) and depth information; take out about 50 grams from each clay sample, put it in an aluminum box of known weight, and use an electronic balance with an accuracy of 0.01 grams to weigh the total weight of the wet soil and the aluminum box. The aluminum box was then placed in an oven at 105-110°C and dried for 6-8 hours until the weight was constant. After being taken out, it was cooled to room temperature in a dryer and the total weight of the dry soil and the aluminum box was weighed again. The moisture content of the sample (the percentage of water in the dry soil weight) was calculated based on the mass difference between the two weighings. The coordinates, depth and corresponding moisture content of each sample were organized into a table to form a discrete moisture content data set. The porosity data measured by the ring knife method (reflecting the proportion of the pore volume of the soil to the total volume) and the moisture content data obtained by the drying method were collected. Each data point contains the corresponding coordinates (X, Y, Z) and Parameter values ​​(porosity, moisture content); divide the three-dimensional grid into regular grids according to the spatial range of the clay area (such as length, width, and height boundaries) and the engineering accuracy requirements (such as a resolution of 0.5 m × 0.5 m × 0.5 m). Each grid node corresponds to a spatial position (X, Y, Z) to be interpolated. For each grid node, search for all known porosity and moisture content data points within a certain range around it (such as a sphere with a radius of 5 meters centered on the node). Usually, multiple points with the closest distance (such as 5-10) are selected as neighboring points to ensure that each node has sufficient neighboring data to support the interpolation calculation.

[0074] According to the principle of "the closer the distance, the greater the influence", the straight-line distance between each neighboring point and the node to be interpolated is calculated, and the reciprocal of the distance (or the square of the reciprocal) is used as the weight of the point. The closer the distance, the higher the weight, and vice versa. The porosity or water content value of the neighboring point is multiplied by the corresponding weight, and the sum is divided by the sum of all weights to obtain the estimated 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 a validation set, compare the interpolation results with the measured values, evaluate the interpolation accuracy (such as calculating the average error and root mean square error), and 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 porosity and water content spatial distribution model.

[0075] The grid nodes are organized into a three-dimensional array by row (X axis), column (Y axis), and layer (Z axis). For example, dataset[i][j][k] corresponds to the node with coordinates (Xi, Yj, Zk), and the porosity n and water content w of the node are stored. For irregular grids or scenes that are convenient for external calls, they can be organized into a table with four or five columns, with each row corresponding to a node. The continuous spatial distribution model of porosity and water content is as follows: Export 3D datasets to formats supported by geotechnical analysis software, such as: CSV files: for general data exchange, can be read by Excel or Python. VTK files: for 3D visualization, support rendering of stereo models in software such as ParaView, Blender, etc. INP files: directly imported into finite element software (such as PLAXIS, FLAC³D) as material parameter input. Extract parameter values ​​of all XY nodes at a specific depth (such as Z=5 meters), generate a plane cloud map of porosity and water content, and use color gradients to represent the value size (such as red for high porosity / water content, blue for low). Through transparency and color mapping, show the continuous change of parameters in space, such as shallow layers (Z=0-5 meters) with higher porosity (semi-transparent red), and deep layers (Z=15-20 meters) with lower porosity (opaque blue).

[0076] The cohesion data (reflecting the bonding ability between clay particles, unit kPa) measured by direct shear test were collected, and the porosity and moisture content values ​​of the corresponding position obtained by interpolation of the above three-dimensional model were matched for each cohesion data point to form a four-dimensional data set containing coordinates, cohesion, porosity, and moisture content.

[0077] Spatial correlation analysis is performed on cohesion, porosity, and water content, and parameter differences at different distances are calculated. For example, the average square of the difference in parameter values ​​of all pairs of points h meters apart is calculated to obtain the variogram value, which reflects the spatial distribution law of the parameter as the distance changes with the distance h (e.g., the closer the distance, the smaller the parameter difference and the stronger the correlation).

[0078] By fitting the variogram curve (such as the spherical model), the spatial correlation range (range), base value (maximum difference) and nugget value (measurement error) of the parameters are determined to quantify the spatial dependence of the parameters. For each grid node to be interpolated, the porosity and water content values ​​of the node are obtained using the known three-dimensional distribution model of porosity and water content. All cohesion data points within the range around the node are searched, and the weighted average cohesion value is calculated by the collaborative Kriging algorithm based on the cohesion values ​​of these points and their porosity and water content differences with the nodes to be interpolated. The algorithm considers both spatial distance and parameter correlation (such as high porosity and high water content areas usually have low cohesion), so that adjacent points with similar parameters contribute more to the interpolation results. 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. The cohesion data is mapped into a spatial distribution map through a visualization tool, which can be presented as a three-dimensional isosurface (such as high cohesion areas are red and low cohesion areas are blue).

[0079] The second source point cloud and the second target point cloud are initially registered, and the global rigid transformation matrix including the rotation matrix and the translation vector is calculated, which specifically 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 and SURF. The feature points of the source point cloud are matched with the feature points of the target point cloud to find the corresponding point pairs. According to the matched point pairs, the initial rotation matrix and translation vector are estimated by the least squares method and other methods. The initial rotation matrix and translation vector are iteratively optimized by the iterative closest point (ICP) algorithm and other methods until the convergence conditions are met to obtain the final global rigid transformation matrix.

[0080] In the embodiment of the present invention, these noise points can be effectively removed by denoising, and the quality of point cloud data can be improved; the denoised point cloud data can be segmented by using the regional growing algorithm, and the second source point cloud and the second target point cloud corresponding to the clay area can be accurately extracted, and the point cloud area with similar characteristics can be automatically identified and extracted, which helps to reduce the amount of calculation for subsequent processing and improve processing efficiency. Downsampling the segmented point cloud can reduce the density of point cloud data; by obtaining clay samples at different depths and the corresponding cohesion parameters, the physical and mechanical properties of clay can be fully understood; based on clay samples and the corresponding cohesion parameters, a three-dimensional distribution model of porosity and water content is established, which can more comprehensively understand the distribution of the physical and mechanical properties of clay in space; the spatial distribution map of cohesion is established by using the collaborative Kriging spatial interpolation method, which 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 accurate registration, the utilization rate of point cloud data can be improved.

[0081] In a preferred embodiment of the present invention, a rigid-non-rigid hybrid registration is performed according to the spatial distribution map and the global rigid transformation matrix to obtain a registered point cloud, including: According to the cohesion value of each point, the corresponding cohesion weight value is calculated, which specifically includes: 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 with coordinates, it is determined that a point is located in a high cohesion area (such as cohesion 50kPa) or a low cohesion area (such as cohesion 20kPa); assigning a "deformation constraint weight" to each point according to the cohesion value, the higher the cohesion point, the greater the weight value (that is, the stronger the constraint on deformation). For example, the weight corresponding to the maximum cohesion value is set to 1 (the strongest constraint), the weight corresponding to the minimum value is set to 0.1 (the weakest constraint), and the intermediate value is assigned by proportional interpolation. 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.

[0082] For each source point after rigid transformation, find the corresponding point in the target point cloud and calculate the current deformation error, which specifically includes: 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, so that it is close to the target point cloud in terms of macroscopic posture, for example, correcting the overall inclination angle or horizontal displacement of the stratum, so that the overall shape of the two point clouds is initially aligned; for each source point after rigid transformation (called "current source point"), search the point with the closest spatial position in the target point cloud as the "corresponding target point". Usually, a spatial neighbor search algorithm (such as KD tree search) is used to find the point with the closest Euclidean distance (the error is within the preset tolerance, such as 5cm), ensuring 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. This difference is the "deformation error", which reflects the local misalignment problem that still exists after rigid transformation. For example, after rigid transformation, a point still has a 2cm deviation in the X direction, no deviation in the Y direction, and a 1cm deviation in the Z direction. The total deformation error is the composite value of these three components.

[0083] According to the cohesion weight value of the source point, the deformation allowable range is adjusted, and the control point position is continuously adjusted by the gradient descent method until the convergence condition is reached to obtain the optimized deformation field, which specifically includes: dynamically setting the "deformation allowable range" according to the cohesion weight value of each source point. Points with high weight (such as cohesion weight 0.9) have low tolerance for deformation errors and only allow small adjustments (such as ±1mm); points with low weight (such as cohesion weight 0.3) have high tolerance and allow larger deformations (such as ±5mm). For example, in the old clay area with high cohesion, any significant local deformation is regarded as abnormal and needs to be strictly corrected; while in the silty clay area with low cohesion, plastic flow deformation within a certain range is allowed; a number of "control points" are uniformly selected in the point cloud (such as selecting a point every 1 meter), and these control points serve as the adjustment hub of the deformation field. Initially, the control point position is the coordinate after rigid transformation. Iterative optimization process: error calculation, for each control point, calculate its deformation error with the corresponding target point, and weight the error according to the weight value (the error contribution of high-weight points is greater). Gradient descent adjustment: adjust the control point position in the opposite direction of the error gradient through the gradient descent method, so that the weighted total error is gradually reduced. Each time the adjustment is made, the control points in the high-weight area move a small amount (to avoid excessive deformation), and the control points in the low-weight area can move a large amount (to adapt to plastic deformation). Convergence judgment: when the total error change of two consecutive iterations is less than the preset threshold (such as 0.5mm), or the maximum number of iterations (such as 100) is reached, stop the adjustment and obtain the optimized deformation field (that is, the final displacement adjustment of each control point).

[0084] Apply the optimized deformation field to the rigidly transformed point cloud to obtain the registered point cloud, which includes: mapping the optimized deformation field (including the displacement adjustment of all control points) to all points of the entire point cloud through interpolation or smoothing algorithm. For example, using thin plate spline interpolation or radial basis function interpolation, the displacement adjustment of non-control points is calculated according to the displacement of the control points to ensure that the deformation field is continuous and conforms to the plastic deformation law of clay (such as the displacement difference of adjacent points changes reasonably with the difference in cohesion weight).

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

[0086] In an embodiment of the present invention, a plastic deformation gradient adjustment function is constructed by combining the spatial distribution map and the global rigid transformation matrix, which can more accurately describe the deformation relationship between point clouds and help to more accurately adjust the position and posture of point clouds during the registration process, thereby improving the accuracy of registration. The rigid-non-rigid hybrid registration method can handle both global rotation and translation transformations (rigid registration) and local deformation and distortion (non-rigid registration), and can more comprehensively consider the deformation between point clouds, thereby improving the accuracy and robustness of registration. According to the cohesion value of each point, the corresponding cohesion weight value is calculated, and the cohesion information can be integrated into the registration process, which can reflect the cohesion difference between point clouds, and help to adaptively adjust the deformation allowable range during the registration process and enhance the adaptability of the model. For each source point after rigid transformation, the corresponding point is found in the target point cloud, and the current deformation error is calculated, so that the quality of registration can be evaluated in real time, and the efficiency of registration can be improved. According to the cohesion weight value of the source point, the deformation allowable range is adjusted, and the degree and direction of deformation can be adaptively controlled, avoiding unnecessary waste of computing resources and improving the efficiency of registration. By continuously adjusting the positions of control points through the gradient descent method, the optimal solution can be gradually approached, and the robustness of the registration 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.

[0087] In a preferred embodiment of the present invention, based on the registered point cloud, calculating the corrected displacement field sequence includes: 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 including: from the registered point cloud data, select point clouds at two different time points (such as time t0 and t1), respectively defined as the reference point cloud (initial state, such as before construction) and the target point cloud (current state, such as after construction), ensure 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 technology (such as KD tree and octree) for fast search, calculate the Euclidean distance and select the point with the smallest distance (the error must be within the preset tolerance, such as 5cm, to ensure point cloud density matching); Example: The target point Pt (10, 20, 5) finds the nearest neighbor point Pb (10.01, 19.98, 5.02) in the reference point cloud and regards it as a deformed point at the same position; For each pair of matching points (Pb, Pt), calculate the three-dimensional displacement vector (△x, △y, △z) = (xt-xb, yt-yb, zt-zb), which reflects the spatial displacement of the target point relative to the reference point (unit: meter), and summarize the displacement vectors of all points to form a data set containing the coordinates of each point and the corresponding displacement.

[0088] 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 including: matching the coordinates (x, y, z) of each point in the initial displacement field with the spatial distribution map of cohesion to obtain the cohesion value c (unit: kPa) corresponding to the point; if the point is located on the grid node of the map, the cohesion value is directly read; if the point is located in the grid unit, 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; example: the coordinates of a point in the displacement field are (12.3, 18.7, 5.5), and the cohesion obtained from the map through interpolation is 45kPa, indicating that the point is located in the medium cohesion clay area.

[0089] According to the magnitude of cohesion, the displacement vector in the initial displacement field is scaled and adjusted to obtain the adjusted displacement vector, including: cohesion reflects the bonding strength between clay particles. High cohesion areas (such as old clay, c>50kPa) have strong anti-deformation ability. The measured displacement may be overestimated due to registration error, and the displacement vector needs to be reduced; low cohesion areas (such as silty clay, c<30kPa) have weak anti-deformation ability and allow larger displacement. The displacement vector can be appropriately enlarged. The scaling factor, s, is set according to the magnitude of cohesion. For example: When c≥50kPa, s=0.8 (reduced by 20% displacement); when 30kPa<c<50kPa, s=1.0 (no adjustment); when c≤30kPa, s=1.2 (enlarged by 20% displacement); for each displacement vector, multiply it by the corresponding scaling factor s to obtain the adjusted displacement vector. Example: The initial displacement of a point is (0.1m, 0, 0), and the cohesion is 25kPa (low cohesion area). The adjusted displacement is (0.1×1.2, 0, 0)=(0.12m, 0, 0), which reflects that the plastic flow deformation of soft clay is reasonably enlarged.

[0090] The adjusted displacement vectors of all points are combined to form a corrected displacement field, the point cloud data at each time point is processed to obtain the corrected displacement field at each time point, and the corrected displacement fields at different time points are arranged in chronological order to generate a corrected displacement field sequence, which specifically includes: integrating the adjusted displacement vectors of all points with the coordinate information to form a corrected displacement field, the data structure is consistent with the initial displacement field, but the displacement vector has been corrected according to the cohesion; For the point cloud data at each time point (such as t0, t1, t2, ..., tn), repeat the above steps (the reference point cloud is updated to the registered point cloud at the previous time point in sequence), and 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 the geological conditions change significantly); the corrected displacement fields of all time points are arranged 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 area, for example: calculate the displacement rate at adjacent time points and identify the deformation acceleration area (such as the displacement rate of the soft soil layer with low cohesion is significantly higher than that of the hard soil layer).

[0091] In the embodiment of the present invention, the deformation of the clay area at different time points can be captured, and the corresponding relationship between the target point cloud and the reference point cloud can be established. This corresponding relationship can ensure that the calculation of the displacement vector is based on the correct point pair, improve the accuracy of the displacement field calculation, and make the displacement vector of each point have its corresponding cohesion value. This matching can consider the influence of cohesion on the displacement vector, adjust the direction and size of the displacement vector according to the spatial distribution of cohesion, so that the displacement field is more consistent with the physical properties of the actual clay area, and improve the accuracy of the displacement field calculation. By considering the constraint effect of cohesion on the displacement vector, the displacement field can be made more consistent with the deformation law of the actual clay area, enhance the adaptability of the model to the physical properties 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 fully reflect the deformation of the clay area at different time points.

[0092] In a preferred embodiment of the present invention, the displacement correction factor and the corrected displacement field sequence are analyzed in the time-frequency domain to extract the characteristic frequency and damping ratio to construct a time-space correlation matrix and generate a displacement data sequence containing the time-space variation characteristics of the geological parameters, including: The displacement correction factor and the corrected displacement field sequence are Fourier transformed to convert them from the time domain to the frequency domain to obtain the representation of the displacement correction factor and the displacement field in the frequency domain, specifically including: obtaining the displacement correction factors of all time points calculated above (each point corresponds to a correction coefficient that varies 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 coordinate + time + displacement component); Fourier transform is performed on the displacement correction factor time series of each spatial point (such as the correction factor of point A from t1 to t100) and the displacement field component time series (such as the displacement value of Δx of point A from t1 to t100), respectively. Each time series is decomposed into a 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. A frequency-amplitude relationship curve in the frequency domain is obtained for each point, 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 the frequency component).

[0093] Analyze the data in the frequency domain, identify the frequency components and the corresponding amplitudes, and in the frequency domain analysis, identify the displacement correction factor and the characteristic frequency of the displacement field sequence; calculate the damping ratio based on the characteristic frequency and the corresponding amplitude, including: Analyze the amplitude spectrum of the frequency domain data and identify the frequency components whose amplitudes are significantly greater than the noise level (i.e., "characteristic frequencies"). For example, find the frequency corresponding to the peak in the amplitude spectrum. This frequency represents the main fluctuation period of deformation in the clay region (e.g., high-frequency components correspond to short-term vibrations, and low-frequency components correspond to long-term settlements). The characteristic frequency of high-cohesion areas may be biased towards low frequencies (slow deformation, such as long-term settlement of foundations), and low-cohesion soft soil layers 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 attenuation of the amplitude with frequency in the frequency domain. The specific steps are: Observe the amplitude attenuation trend near the characteristic frequency. The amplitude in areas with high damping ratio (such as soft clay) decreases faster with increasing frequency, while the amplitude in areas with low damping ratio (such as hard clay) decays slowly. For each characteristic frequency, compare the measured amplitude with the theoretical amplitude of undamped vibration. Estimate the damping ratio (dimensionless, usually between 0-1) corresponding to the frequency through the slope of the attenuation curve or the amplitude ratio of adjacent peaks. The higher the damping ratio, the stronger the energy absorption capacity of the clay and the faster the deformation fluctuation decays (for example, the damping ratio of silty clay can reach 0.2-0.3, while the damping ratio of old clay is only 0.05-0.1).

[0094] The extracted characteristic frequencies and damping ratios are combined with the spatial and temporal information of the displacement correction factors and displacement field sequences to construct a spatiotemporal correlation matrix, which includes: extracting the following information for each spatial point (x, y, z) and each time point t: spatial information: coordinates, cohesion value, porosity and other geological parameters; time information: timestamp, sampling interval (such as daily / weekly monitoring); frequency domain features: characteristic frequencies (such as f1, f2), corresponding damping ratios (ξ1, ξ2), and amplitudes of each frequency component. Matrix structure design: Row dimension, spatial points (arranged by grid nodes or discrete points); column dimension, time points + frequency domain features (such as time t1-tn, characteristic frequency f1-fm, damping ratio ξ1-ξm).

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

[0096] Based on the spatiotemporal correlation matrix, the original displacement field sequence is corrected to generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters, including: adjusting the displacement field sequence at each time point according to the frequency domain characteristics in the spatiotemporal correlation matrix: frequency component screening, retaining characteristic frequency components that are 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 areas with high damping ratios (such as the displacement amplitude of soft soil layers at seismic frequencies needs to be multiplied by the damping correction factor). Position correction, considering the phase difference of characteristic frequencies at different spatial points (such as the displacement phase lag caused by stratum tilt), adjusting the phase parameters of the displacement vector, so that the deformation field is more consistent with geological continuity in time and space (such as the displacement phase difference of adjacent points at the same frequency does not exceed the preset threshold).

[0097] The spatial distribution of geological parameters such as cohesion and porosity is used as weight to weight the displacement data sequence: the displacement sequence in the high cohesion area is dominated by low-frequency components and the displacement amplitude is small; the low cohesion area allows more high-frequency components (such as deformation fluctuations caused by seepage) and the amplitude is large. For example, in the hard clay area with a cohesion of 50kPa, the high-frequency (>1Hz) components in the corrected displacement sequence are suppressed, and the low-frequency (<0.1Hz) settlement signal is retained and amplified; while in the soft soil area with a cohesion of 20kPa, the medium-frequency (0.5-1Hz) vibration component is retained, reflecting its characteristics that are easily affected by dynamic loads. The corrected displacement vectors of all spatial points are integrated and arranged in chronological order to form a displacement data sequence containing spatiotemporal variation characteristics. The data at each time point includes: the coordinates of each point, geological parameters such as cohesion and porosity; the corrected three-dimensional displacement vector (considering the adjustment of characteristic frequency and damping ratio); the amplitude and phase information of each frequency component.

[0098] In the embodiment of the present invention, the displacement correction factor and the displacement field sequence are converted from the time domain to the frequency domain by Fourier transform. In the frequency domain, the amplitude distribution of the displacement field on different frequency components can be clearly seen; in the frequency domain analysis, the characteristic frequency of the displacement correction factor and the displacement field sequence is 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 of the displacement field during the vibration process. The extracted characteristic frequency and damping ratio are combined with the spatial and temporal information of the displacement correction factor and the displacement field sequence to construct a spatiotemporal correlation matrix. The spatiotemporal correlation matrix can reflect the correlation between the displacement correction factor and the displacement field sequence in space and time, and is helpful to reveal the influence of the spatiotemporal variation characteristics of geological parameters on the displacement field; based on the spatiotemporal correlation matrix, the original displacement field sequence is corrected to generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters, and the displacement field can be dynamically adjusted according to the changes in geological parameters, thereby improving the accuracy of the displacement data sequence. The displacement data sequence containing the spatiotemporal variation characteristics of geological parameters can reflect the impact of changes in geological parameters on the displacement field. By comparing and analyzing the displacement data sequences under different geological parameter conditions, the deformation under different geological conditions can be predicted.

[0099] Among them, after the above step 3, the generated displacement data sequence containing the temporal and spatial variation characteristics of geological parameters can also be input into the monitoring system to realize real-time dynamic deformation monitoring of rainwater resource interception structures, and by comparing the displacement data at different time nodes, the subtle deformation trends of the structures under complex geological conditions can be accurately captured.

Claims

1. An intelligent monitoring method for displacement of rainwater resource interception and storage structures, characterized in that: The method comprises: Based on the initial displacement field constraints obtained by point cloud registration, a three-dimensional spatial distribution model is constructed by coupling multi-source geological data through particle swarm optimization; the sandy soil area is discretized into dynamic subdomain grids, and parameter values ​​are assigned to each subdomain based on the three-dimensional spatial distribution model to form a parameterized grid; based on the parameterized grid, an elastic modulus correction model is established; based on the elastic modulus correction model, the deformation field of the particle swarm is iterated to obtain the adjusted deformation field; based on the adjusted deformation field, the displacement correction factor is calculated; Obtain the second source point cloud and the second target point cloud of the clay area and establish a spatial distribution map; 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 map and the global rigid transformation matrix to obtain the registered point cloud; calculate the corrected displacement field sequence based on the registered point cloud; The displacement correction factors and the corrected displacement field sequences are analyzed in the time-frequency domain to extract characteristic frequencies and damping ratios to construct a spatiotemporal correlation matrix and generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters.

2. According to claim 1, a method for intelligently monitoring the displacement of a rainwater resource interception and storage structure is characterized in that: Based on the initial displacement field constraints obtained by point cloud registration, a three-dimensional spatial distribution model is constructed by coupling multi-source geological data through particle swarm optimization, including: Obtain source point cloud and target point cloud of sandy soil area, and obtain discrete geological exploration data, including permeability coefficient, porosity and 3D coordinates of the measuring points; Align the denoised source point cloud and target point cloud and calculate the initial displacement field; The geological exploration data of each measuring point is converted into a particle, and the particle attributes include permeability coefficient, porosity and spatial coordinates; a particle group is generated in three-dimensional space, and the spatial correlation of geological parameters is simulated through virtual "force"; Calculate the evaluation value of the geological parameters of the current particle position, update the particle speed and position, and stop the optimization when the particle group converges to obtain the converged particle group; Based on the converged particle swarm, a continuous three-dimensional spatial distribution model of permeability and porosity is generated.

3. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 2 is characterized in that: Calculate the estimated values ​​of the geological parameters at the current particle position, including: Calculate the theoretical permeability coefficient at the particle location based on the current particle porosity and the effective particle size of the soil particles; Compare the current particle measured permeability coefficient with the theoretical permeability coefficient, calculate the absolute value of the relative deviation, and obtain the physical constraint value; Taking the current particle as the center, determine the set of neighboring particles around it according to the preset neighborhood search radius, and calculate the standard deviation of permeability coefficient and porosity within the neighborhood search radius; calculate the spatial constraint value according to the standard deviation of permeability coefficient and porosity; The physical constraint value and the spatial constraint value are weighted and summed to obtain the evaluation value.

4. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 3 is characterized in that: The sandy soil area is discretized into dynamic subdomain grids, and parameter values ​​are assigned 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 grid data and the geometric information of the soil layer interface, the subdomains are divided step by step on the two-dimensional plane to obtain the dynamic subdomain grid; According to the permeability and porosity values ​​of the dynamic subdomain grid and the corresponding grid cells, each subdomain is analyzed to calculate the average and maximum values ​​of the permeability, the average and minimum values ​​of the porosity, and the coefficient of variation of the permeability and porosity within the subdomain; The average value, maximum value, minimum value and coefficient of variation of the subdomain are combined into a parameter vector as the geological attribute of the subdomain; Normalizing the parameter vector to obtain a normalized parameter vector; The normalized parameter vector is used as the geological attribute parameter and appended to the corresponding grid cell to form a parameterized grid. Each grid cell contains the vertex coordinates and the corresponding geological attributes.

5. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 4 is characterized in that: According to the parameterized grid, an elastic modulus correction model is established; According to the elastic modulus correction model, the deformation field of the particle group is iterated to obtain the adjusted deformation field, including: Determine the joint distribution constraint relationship between porosity and permeability; Obtain the measured data of elastic modulus and porosity of samples of different lithologies, analyze their corresponding relationship, and establish the lithological function relationship between elastic modulus and porosity for different lithologies; According to the joint distribution constraint relationship and lithology function relationship, the elastic modulus correction model of joint constraint of permeability coefficient and porosity is calculated; According to the elastic modulus correction model, the permeability coefficient and porosity data of the parameterized grid, the static elastic modulus value of each grid unit is calculated to form an initial elastic modulus field; The elastic modulus values ​​of all grid cells are combined into a three-dimensional spatial distribution model, and the simulated deformation field is calculated; Compare the simulated deformation field with the actual observed data to obtain difference analysis results; According to the results of the difference analysis, the deformation field parameters are adjusted to obtain the adjusted deformation field; according to the adjusted deformation field, the pore pressure field changes are inverted, and the elastic modulus field is updated through the pore pressure changes. When the deformation field adjustment amount is less than the preset threshold, the current round of optimization is stopped. When the difference in the deformation field between two consecutive iterations is lower than the preset threshold, the final deformation field and elastic modulus field are output.

6. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 5 is characterized in that: Calculate the displacement correction factor based on the adjusted deformation field, including: According to the geological stratification information of the sandy soil area, the deformation field data is divided into several horizontal layers; For each layer, the displacement vectors of all points in the layer are calculated, and the characteristic displacement of the layer is calculated. The characteristic displacement includes the maximum displacement, the minimum displacement and the average displacement; Obtain the sand state parameters of the sandy soil area, and spatially match the sand state parameters with the deformation field data so that the displacement vector of each point has its corresponding sand state parameter; Using characteristic displacement and sand state parameters, calculate the displacement correction coefficient of each layer; The displacement correction factor is calculated based on the adjusted deformation field and the displacement correction coefficient.

7. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 1 is 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: Original point cloud data is collected in the clay area, and denoising is performed on the original point cloud to obtain denoised point cloud data; based on the denoised point cloud data, a second source point cloud and a second target point cloud corresponding to the clay area are segmented by a regional growing algorithm, and the segmented point cloud is downsampled to obtain a second source point cloud and a second target point cloud; Obtain clay samples at different depths and the corresponding cohesion parameters; establish a three-dimensional distribution model of porosity and water content based on the clay samples and the corresponding cohesion parameters; establish a cohesion spatial distribution map using the collaborative Kriging spatial interpolation method based on the three-dimensional distribution model of porosity and water content; The second source point cloud and the second target point cloud are initially registered, and a global rigid transformation matrix including a rotation matrix and a translation vector is calculated.

8. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 7 is characterized in that: According to the spatial distribution map and the global rigid transformation matrix, a rigid-non-rigid hybrid registration is performed to obtain the registered point cloud, including: According to the cohesion value of each point, the corresponding cohesion weight value is calculated; For each source point after rigid transformation, find the corresponding point in the target point cloud and calculate the current deformation error; According to the cohesion weight value of the source point, the deformation allowable range is adjusted, and the control point position is continuously adjusted by the gradient descent method until the convergence condition is reached to obtain the optimized deformation field; The optimized deformation field is applied to the rigidly transformed point cloud to obtain the registered point cloud.

9. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 8 is characterized in that: Based on the registered point cloud, the corrected displacement field sequence is calculated, including: 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; 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; According to the magnitude of the cohesive force, the displacement vector in the initial displacement field is scaled and adjusted to obtain an adjusted displacement vector; The adjusted displacement vectors of all points are combined to form a corrected displacement field, and the point cloud data at each time point is processed to obtain the corrected displacement field at each time point; The corrected displacement fields at different time points are arranged in time sequence to generate a corrected displacement field sequence.

10. The intelligent monitoring method for displacement of rainwater resource interception and storage structures according to claim 9 is characterized in that: The displacement correction factor and the corrected displacement field sequence are analyzed in the time-frequency domain to extract the characteristic frequency and damping ratio to construct a spatiotemporal correlation matrix and generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters, including: Performing Fourier transformation on the displacement correction factor and the corrected displacement field sequence, converting them from the time domain to the frequency domain, so as to obtain the representation 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 in the frequency domain analysis, identify the displacement correction factor and the characteristic frequency of the displacement field sequence; calculate the damping ratio based on the characteristic frequency and the corresponding amplitude; The extracted characteristic frequencies and damping ratios are combined with the displacement correction factors and the spatial and temporal information of the displacement field sequence to construct a spatiotemporal correlation matrix; Based on the spatiotemporal correlation matrix, the original displacement field sequence is modified to generate a displacement data sequence containing the spatiotemporal variation characteristics of geological parameters.

Images