A foundation pit deformation monitoring method and system

Abstract

The application provides a foundation pit deformation monitoring method and system, which comprises the following steps: acquiring a three-dimensional coordinate time series of GNSS monitoring points around a foundation pit, calculating a three-dimensional velocity field, estimating a velocity gradient tensor, extracting a principal strain rate based on a symmetric strain rate component of the velocity gradient tensor, establishing a nonlinear mapping of the principal strain rate and a neighborhood search influence factor, adjusting a density clustering core point discrimination condition, generating multiple candidate regional division schemes, constructing a joint penalty function containing a residual velocity discontinuity norm and a model complexity penalty term, selecting an optimal regional division scheme, combining observation errors and model prediction residuals, and calculating local translation parameters and rotation parameters of each region by using an iterative reweighted least squares method, and then performing topological optimization according to the differences between the motion parameters of the regions, merging adjacent regions that meet the motion consistency condition, and reattributing boundary monitoring points until the regional division is stable, and finally outputting a foundation pit deformation field model.

Description

Technical Field

[0001] This application belongs to the field of foundation pit deformation monitoring, and in particular relates to a method and system for foundation pit deformation monitoring. Background Technology

[0002] Foundation pit engineering is a crucial component of urban underground space development and large-scale building construction. Its stability directly impacts construction progress, on-site operational safety, and the safe operation of surrounding existing buildings, roads, pipelines, and underground facilities. During foundation pit excavation and support construction, the surrounding strata are prone to various deformation phenomena, such as settlement, horizontal displacement, and localized shear deformation, due to the combined influence of geological conditions, groundwater level changes, the stress state of the support structure, and construction disturbances. To promptly grasp the evolution trend of the foundation pit's safety status and reduce risks such as collapse, cracking of surrounding buildings, and damage to underground pipelines, continuous and accurate monitoring of foundation pit deformation is of paramount importance. In recent years, Global Navigation Satellite Systems (GNSS) have gained increasing traction due to their advantages, including all-weather observation, continuous acquisition of three-dimensional coordinate time series data, and suitability for long-term automated monitoring. Deploying a GNSS monitoring network around the foundation pit allows for the continuous acquisition of spatial location change data from multiple monitoring points, providing a data foundation for foundation pit deformation analysis.

[0003] However, due to the significant spatial differences in ground stress release, support structure constraints, and groundwater environment changes during foundation pit excavation, the actual deformation of the foundation pit is usually not a rigid body motion, but rather exhibits spatial heterogeneity, local discontinuity, and staged evolution characteristics. Existing GNSS deformation monitoring data processing methods often rely on manually set empirical parameters or use fixed global motion models for region division and parameter calculation when constructing deformation field models, making it difficult to adequately adapt to the different motion states between different blocks. When discontinuous sliding, differential settlement, or slight torsion occur in local areas of the foundation pit, traditional processing methods are prone to problems such as unclear deformation boundary identification, large deviations in local motion parameter estimation, and oversimplification or overfitting of the model, thus affecting the ability of the deformation field model to express the actual local motion characteristics. Therefore, how to extract strain characteristics that reflect local deformation differences from the three-dimensional velocity field formed by GNSS monitoring points, reduce the influence of fixed empirical parameters on the region division results, and achieve adaptive identification, optimized division, and robust calculation of motion parameters for local deformation regions is a technical problem that urgently needs to be solved in the field of foundation pit safety monitoring data processing. Summary of the Invention

[0004] To address the problems in existing technologies that rely on fixed empirical parameters and global motion models, making it difficult to adapt to local variations in foundation pit deformation, resulting in unclear regional boundaries, large deviations in local motion parameter estimation, and inaccurate expression of deformation field models, this invention proposes a foundation pit deformation monitoring method and system.

[0005] In a first aspect, the present invention proposes a method for monitoring the deformation of a foundation pit, comprising the following steps:

[0006] The three-dimensional coordinate time series of GNSS monitoring points around the foundation pit is obtained, and the three-dimensional velocity field of each monitoring point is calculated. The velocity gradient tensor is estimated based on the three-dimensional velocity field, and the principal strain rate representing the degree of deformation is extracted according to the symmetric strain rate component of the velocity gradient tensor. The principal strain rate corresponding to the triangular unit is mapped to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point. Through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor, the core point discrimination condition of the density clustering algorithm is adjusted. The updated core point discrimination condition is applied to the three-dimensional velocity field for clustering processing to obtain multiple candidate region division schemes.

[0007] For each candidate region partitioning scheme, a joint penalty function is constructed. All candidate region partitioning schemes are traversed to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme. For each region of the optimal region partitioning scheme, the local translation parameters and local rotation parameters of each region are calculated using the iterative reweighted least squares method based on the mean square error of the observation values ​​of each monitoring point and the model prediction residual.

[0008] Based on the local translation parameters and local rotation parameters of each region, a topology optimization iterative operation is performed. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency condition, the adjacent regions are merged. The boundary monitoring points with smaller fitting residuals for adjacent regions are reassigned. The region merging and boundary monitoring point reassignment are alternately executed until the region division is stable, and the foundation pit deformation field model is output.

[0009] Optionally, the step of obtaining the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit and calculating the three-dimensional velocity field of each monitoring point includes:

[0010] The three-dimensional coordinate time series is denoised and gross error removed to obtain a smoothed three-dimensional coordinate time series.

[0011] The smoothed three-dimensional coordinate time series was linearly fitted using the sliding window method, and the deformation rate within each time window was calculated.

[0012] The deformation rate within each time window is time-weighted and averaged to obtain the three-dimensional velocity of each monitoring point in the north, east and vertical directions, and combined to obtain the three-dimensional velocity field.

[0013] Optionally, the step of estimating the velocity gradient tensor based on the three-dimensional velocity field, extracting the principal strain rate representing the degree of deformation based on the symmetric strain rate component of the velocity gradient tensor, and mapping the principal strain rate corresponding to the triangular element to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point includes:

[0014] A Delaunay triangulation was constructed with all GNSS monitoring points as vertices, and triangular units that crossed the cavities inside the pit were removed to obtain the continuous topology around the pit.

[0015] For each triangular element of the Delaunay triangulation, a system of partial differential equations is established based on the planar coordinates of the vertices and the corresponding horizontal velocity components, and the two-dimensional velocity gradient tensor of each element is obtained by solving the equations.

[0016] The symmetric part of the two-dimensional velocity gradient tensor and its transpose matrix is ​​obtained to obtain the two-dimensional strain rate tensor. The eigenvalues ​​of the two-dimensional strain rate tensor are calculated, and the eigenvalue with the largest absolute value is extracted as the principal strain rate of the triangular element.

[0017] The principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are statistically fused to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point.

[0018] Optionally, the nonlinear mapping from the principal strain rate to the neighborhood search influence factor is used to adjust the core point discrimination condition of the density clustering algorithm, and the updated core point discrimination condition is applied to cluster the three-dimensional velocity field to obtain multiple candidate region partitioning schemes, including:

[0019] An exponential decay function with the dimensionless principal strain rate as the independent variable is constructed, and a neighborhood search influence factor between 0 and 1 is calculated.

[0020] The minimum number of contained points preset in the density clustering algorithm is multiplied by the neighborhood search influence factor to obtain the updated core point discrimination threshold. The larger the principal strain rate, the smaller the neighborhood search influence factor, and the lower the corresponding core point discrimination threshold.

[0021] Given a neighborhood radius, the three-dimensional velocity field is subjected to multiple density clustering scans using an updated core point discrimination threshold to obtain multiple candidate region partitioning schemes.

[0022] Optionally, the step of constructing a joint penalty function for each candidate region partitioning scheme, and traversing all candidate region partitioning schemes to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, includes:

[0023] The validity of the sub-regions in each candidate region partitioning scheme is determined, and candidate region partitioning schemes that do not meet the preset minimum number of solution points or preset rank conditions are eliminated.

[0024] The residual velocity discontinuity norm of all monitoring points on the boundary of adjacent regions in each candidate region division scheme is calculated as a data fitting term;

[0025] A model complexity penalty term is constructed based on the number of regions. The boundary principal strain rate statistics are calculated according to the principal strain rate of the monitoring points corresponding to the boundary monitoring points. The boundary principal strain rate statistics are then dimensionless to generate dimensionless reduction weights. These dimensionless reduction weights are then used to reduce the model complexity penalty term.

[0026] The data fitting term after being dimensionless by the benchmark velocity norm is added to the reduced dimensionless model complexity penalty term to form the joint penalty function. All candidate region partitioning schemes are traversed, and the scheme with the smallest value of the joint penalty function is selected as the optimal region partitioning scheme. The number of regions at this time is the optimal number of local motion models.

[0027] Optionally, the step of calculating the local translation and local rotation parameters of each region using the iterative reweighted least squares method based on the mean square error of the observed values ​​at each monitoring point and the model prediction residual includes:

[0028] An initial diagonal weight matrix is ​​constructed using the reciprocal of the squared error of the observation values ​​at each monitoring point, and the initial local translation parameters and initial local rotation parameters are initially calculated.

[0029] The model-predicted velocity is calculated based on the initial local translation parameters and the initial local rotation parameters, and the vector residual between the model-predicted velocity and the actual observed velocity is obtained.

[0030] The standardized residuals are calculated based on the vector residuals and the unit weight mean square error of the overall adjustment, and a robust weight function is constructed based on the standardized residuals to obtain a new weight matrix.

[0031] The new weight matrix is ​​used to replace the original weight matrix for iterative solution until the change in local translation parameters between two consecutive solutions is less than the preset translation convergence threshold, and the change in local rotation parameters between two consecutive solutions is less than the preset rotation convergence threshold. The local translation parameters and local rotation parameters are then output.

[0032] Optionally, the step of performing topology optimization iteration based on the local translation parameters and local rotation parameters of each region, merging adjacent regions when the differences in translation and rotation parameters between adjacent regions meet a preset motion consistency condition, and reassigning the boundary monitoring points with smaller fitting residuals to adjacent regions, alternately performing region merging and boundary monitoring point reassignment until the region division is stable, includes:

[0033] Calculate the difference in local translation vectors and the difference in local rotation vectors between two adjacent regions. Based on the average enclosing radius of the two adjacent regions, convert the difference in local rotation vectors into an equivalent linear velocity difference. Calculate the composite difference norm based on the difference in local translation vectors and the equivalent linear velocity difference.

[0034] If the composite difference norm is lower than the preset composite consistency threshold, then the two adjacent regions are merged, and the local translation parameters and local rotation parameters of the merged region are recalculated.

[0035] Calculate the fitting residual of the boundary monitoring point to the local motion model of the original region to which the monitoring point belongs, and the fitting residual of the local motion model of all adjacent regions. If the fitting residual of a certain adjacent region is less than the fitting residual of the original region, and the difference between the two satisfies the preset residual reduction condition, then the monitoring point is reclassified to the adjacent region.

[0036] The process alternately executes the region merging and boundary monitoring point reassignment until no adjacent region meets the preset motion consistency condition and no monitoring point undergoes a change in assignment. Then, the iteration terminates and the foundation pit deformation field model is output.

[0037] On the other hand, the present invention also proposes a foundation pit deformation monitoring system, comprising the following modules:

[0038] The processing module is used to acquire the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit, calculate the three-dimensional velocity field of each monitoring point, estimate the velocity gradient tensor based on the three-dimensional velocity field, extract the principal strain rate representing the degree of deformation according to the symmetric strain rate component of the velocity gradient tensor, map the principal strain rate corresponding to the triangular unit to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point, adjust the core point discrimination condition of the density clustering algorithm through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor, apply the updated core point discrimination condition to cluster the three-dimensional velocity field, and obtain multiple candidate region division schemes.

[0039] The module is used to construct a joint penalty function for each candidate region partitioning scheme, traverse all candidate region partitioning schemes and select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, and for each region of the optimal region partitioning scheme, calculate the local translation parameters and local rotation parameters of each region based on the mean square error of the observation values ​​of each monitoring point and the model prediction residual using the iterative reweighted least squares method.

[0040] The output module is used to perform topology optimization iteration based on the local translation parameters and local rotation parameters of each region. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency conditions, the adjacent regions are merged. The boundary monitoring points with smaller fitting residuals to the adjacent regions are reassigned. The region merging and boundary monitoring point reassignment are alternately executed until the region division is stable, and the foundation pit deformation field model is output.

[0041] Preferably, the step of acquiring the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit and calculating the three-dimensional velocity field of each monitoring point includes:

[0042] The three-dimensional coordinate time series is denoised and gross error removed to obtain a smoothed three-dimensional coordinate time series.

[0043] The smoothed three-dimensional coordinate time series was linearly fitted using the sliding window method, and the deformation rate within each time window was calculated.

[0044] The deformation rate within each time window is time-weighted and averaged to obtain the three-dimensional velocity of each monitoring point in the north, east and vertical directions, and combined to obtain the three-dimensional velocity field.

[0045] Preferably, the step of estimating the velocity gradient tensor based on the three-dimensional velocity field, extracting the principal strain rate representing the degree of deformation based on the symmetric strain rate component of the velocity gradient tensor, and mapping the principal strain rate corresponding to the triangular element to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point includes:

[0046] A Delaunay triangulation was constructed with all GNSS monitoring points as vertices, and triangular units that crossed the cavities inside the pit were removed to obtain the continuous topology around the pit.

[0047] For each triangular element of the Delaunay triangulation, a system of partial differential equations is established based on the planar coordinates of the vertices and the corresponding horizontal velocity components, and the two-dimensional velocity gradient tensor of each element is obtained by solving the equations.

[0048] The symmetric part of the two-dimensional velocity gradient tensor and its transpose matrix is ​​obtained to obtain the two-dimensional strain rate tensor. The eigenvalues ​​of the two-dimensional strain rate tensor are calculated, and the eigenvalue with the largest absolute value is extracted as the principal strain rate of the triangular element.

[0049] The principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are statistically fused to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point.

[0050] Preferably, the nonlinear mapping from the principal strain rate to the neighborhood search influence factor is used to adjust the core point discrimination condition of the density clustering algorithm, and the updated core point discrimination condition is applied to cluster the three-dimensional velocity field to obtain multiple candidate region partitioning schemes, including:

[0051] An exponential decay function with the dimensionless principal strain rate as the independent variable is constructed, and a neighborhood search influence factor between 0 and 1 is calculated.

[0052] The minimum number of contained points preset in the density clustering algorithm is multiplied by the neighborhood search influence factor to obtain the updated core point discrimination threshold. The larger the principal strain rate, the smaller the neighborhood search influence factor, and the lower the corresponding core point discrimination threshold.

[0053] Given a neighborhood radius, the three-dimensional velocity field is subjected to multiple density clustering scans using an updated core point discrimination threshold to obtain multiple candidate region partitioning schemes.

[0054] Preferably, the step of constructing a joint penalty function for each candidate region partitioning scheme, and traversing all candidate region partitioning schemes to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, includes:

[0055] The validity of the sub-regions in each candidate region partitioning scheme is determined, and candidate region partitioning schemes that do not meet the preset minimum number of solution points or preset rank conditions are eliminated.

[0056] The residual velocity discontinuity norm of all monitoring points on the boundary of adjacent regions in each candidate region division scheme is calculated as a data fitting term;

[0057] A model complexity penalty term is constructed based on the number of regions. The boundary principal strain rate statistics are calculated according to the principal strain rate of the monitoring points corresponding to the boundary monitoring points. The boundary principal strain rate statistics are then dimensionless to generate dimensionless reduction weights. These dimensionless reduction weights are then used to reduce the model complexity penalty term.

[0058] The data fitting term after being dimensionless by the benchmark velocity norm is added to the reduced dimensionless model complexity penalty term to form the joint penalty function. All candidate region partitioning schemes are traversed, and the scheme with the smallest value of the joint penalty function is selected as the optimal region partitioning scheme. The number of regions at this time is the optimal number of local motion models.

[0059] Preferably, the step of calculating the local translation and local rotation parameters of each region using the iterative reweighted least squares method based on the mean square error of the observed values ​​at each monitoring point and the model prediction residual includes:

[0060] An initial diagonal weight matrix is ​​constructed using the reciprocal of the squared error of the observation values ​​at each monitoring point, and the initial local translation parameters and initial local rotation parameters are initially calculated.

[0061] The model-predicted velocity is calculated based on the initial local translation parameters and the initial local rotation parameters, and the vector residual between the model-predicted velocity and the actual observed velocity is obtained.

[0062] The standardized residuals are calculated based on the vector residuals and the post-correction unit weight mean square error of the overall adjustment, and a robust weight function is constructed based on the standardized residuals to obtain a new weight matrix.

[0063] The new weight matrix is ​​used to replace the original weight matrix for iterative solution until the change in local translation parameters between two consecutive solutions is less than the preset translation convergence threshold, and the change in local rotation parameters between two consecutive solutions is less than the preset rotation convergence threshold. The local translation parameters and local rotation parameters are then output.

[0064] Preferably, the step of performing topology optimization iteration based on the local translation parameters and local rotation parameters of each region, merging adjacent regions when the differences in translation and rotation parameters between adjacent regions meet a preset motion consistency condition, and reassigning the boundary monitoring points with smaller fitting residuals to adjacent regions, alternately performing region merging and boundary monitoring point reassignment until the region division is stable, includes:

[0065] Calculate the difference in local translation vectors and the difference in local rotation vectors between two adjacent regions. Based on the average enclosing radius of the two adjacent regions, convert the difference in local rotation vectors into an equivalent linear velocity difference. Calculate the composite difference norm based on the difference in local translation vectors and the equivalent linear velocity difference.

[0066] If the composite difference norm is lower than the preset composite consistency threshold, then the two adjacent regions are merged, and the local translation parameters and local rotation parameters of the merged region are recalculated.

[0067] Calculate the fitting residual of the boundary monitoring point to the local motion model of the original region to which the monitoring point belongs, and the fitting residual of the local motion model of all adjacent regions. If the fitting residual of a certain adjacent region is less than the fitting residual of the original region, and the difference between the two satisfies the preset residual reduction condition, then the monitoring point is reclassified to the adjacent region.

[0068] The process alternately executes the region merging and boundary monitoring point reassignment until no adjacent region meets the preset motion consistency condition and no monitoring point undergoes a change in assignment. Then, the iteration terminates and the foundation pit deformation field model is output.

[0069] This invention solves the three-dimensional velocity field of GNSS monitoring points and extracts the principal strain rates. It then uses nonlinear mapping to adaptively adjust the discrimination conditions for density clustering core points, enabling the region partitioning process to incorporate local deformation intensity and reducing the risk of mispartitioning caused by fixed clustering parameters. By constructing a joint penalty function that considers both the residual velocity discontinuity norm and the model complexity penalty term, a better region partitioning scheme is selected, balancing model fitting effect and partitioning complexity. Furthermore, iterative reweighted least squares is used to solve for local translation and rotation parameters, mitigating the impact of abnormal observation errors on parameter estimation. Iterative topology optimization merges adjacent regions and reassigns boundary monitoring points, reducing redundant partitioning and boundary ambiguity, and improving the spatial continuity and local motion representation capability of the foundation pit deformation field model. Attached Figure Description

[0070] Figure 1 A flowchart of the first embodiment;

[0071] Figure 2 A schematic diagram of the fitting curve between the time-series displacement observations of GNSS monitoring points and the linear deformation trend;

[0072] Figure 3 This is a schematic diagram of the horizontal velocity vector field distribution at monitoring points around the foundation pit. Detailed Implementation

[0073] The exemplary embodiments of this disclosure are described below with reference to the accompanying drawings, including various details of the embodiments to aid understanding, and should be considered merely exemplary. Therefore, those skilled in the art will recognize that various changes and modifications can be made to the embodiments described herein without departing from the scope of this disclosure. Similarly, for clarity and brevity, descriptions of well-known functions and structures are omitted in the following description.

[0074] In the first embodiment, the present invention proposes a method for monitoring foundation pit deformation, such as... Figure 1 As shown, it includes the following steps:

[0075] S1. Obtain the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit, calculate the three-dimensional velocity field of each monitoring point, estimate the velocity gradient tensor based on the three-dimensional velocity field, and extract the principal strain rate representing the degree of deformation according to the symmetric strain rate component of the velocity gradient tensor. Map the principal strain rate corresponding to the triangular unit to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point. Adjust the core point discrimination condition of the density clustering algorithm through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor. Apply the updated core point discrimination condition to cluster the three-dimensional velocity field to obtain multiple candidate region division schemes.

[0076] The observation and navigation files of each monitoring point around the foundation pit were acquired using a GNSS receiver. The single-point positioning and carrier phase differential positioning modules in RTKLIB were then used for calculation to obtain the three-dimensional coordinate time series of each monitoring point in the WGS84 coordinate system, oriented northeast. The fitting curve between the displacement time series observations of the GNSS monitoring points and the linear deformation trend is shown below. Figure 2 As shown in the figure, this diagram illustrates the difference between the original noisy time series and the linear trend extraction effect. Using Hector time series analysis software combined with the maximum likelihood estimation method, steps and outliers in the coordinate time series were removed, and a linear trend term was extracted to obtain the three-dimensional velocity field at each monitoring point.

[0077] The monitoring points are divided into Delaunay triangulations, and a system of spatial partial derivative equations of velocity and coordinates is constructed for each triangular element to solve for the velocity gradient tensor. The symmetric portion of the velocity gradient tensor and its transpose is obtained to yield a two-dimensional strain rate tensor. The eigenvalue with the largest absolute value is extracted using eigenvalue decomposition as the principal strain rate of the triangular element. Subsequently, the principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are weighted and averaged according to the area of ​​the triangular elements to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point; where the larger the area of ​​the triangular element, the greater its corresponding weight.

[0078] A monotonically decreasing nonlinear mapping model was constructed. The principal strain rates of each monitoring point were normalized and input into the model to calculate the corresponding neighborhood search influence factor. This influence factor was then multiplied by the preset minimum containment point threshold in the density clustering algorithm to obtain the updated core point discrimination threshold for each monitoring point. Finally, a parameter sequence containing different basic neighborhood radii and minimum containment points was set. Based on the improved DBSCAN density clustering algorithm, the three-dimensional velocity field was subjected to multiple clustering processes, and the connected region results obtained each time were saved as candidate region partitioning schemes.

[0079] In an optional embodiment, the step of acquiring the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit and calculating the three-dimensional velocity field of each monitoring point includes:

[0080] The three-dimensional coordinate time series is denoised and gross error removed to obtain a smoothed three-dimensional coordinate time series.

[0081] The smoothed three-dimensional coordinate time series was linearly fitted using the sliding window method, and the deformation rate within each time window was calculated.

[0082] The deformation rate within each time window is time-weighted and averaged to obtain the three-dimensional velocity of each monitoring point in the north, east and vertical directions, and combined to obtain the three-dimensional velocity field.

[0083] The Symlet 8 wavelet basis was used to decompose the 3D coordinate time series of GNSS monitoring points around the foundation pit into 3 to 5 levels for noise reduction. Gross errors were identified based on the 3σ criterion, and outliers deviating from the median of the local window by more than 3 times the standard deviation were removed. Cubic spline interpolation was used to fill in missing data, resulting in a continuous and smooth 3D coordinate time series. The sliding window was set to a time span of 30 days with a step size of 1 day. Univariate linear regression models were constructed for the northward, eastward, and vertical time series of each monitoring point within the window. The slope parameter of the regression line was solved using the least squares method, and this slope was extracted as the deformation rate within the current window period. For example, the northward deformation rate was found to be 2.5 mm / yr.

[0084] Weighting coefficients are constructed using the reciprocal of the sum of squared residuals fitted for each time window; windows with smaller residuals have higher weights. The deformation rates of all sliding windows generated within 180 days are then time-weightedly averaged according to these weighting coefficients to obtain a more robust final three-dimensional velocity vector for each GNSS monitoring point, for example, 1.2 mm / yr northward, -0.8 mm / yr eastward, and -3.5 mm / yr vertically. The velocity vectors of all monitoring points in the entire perimeter monitoring network of the foundation pit are merged and arranged according to their spatial relationships to construct a three-dimensional velocity field representing the overall surface movement trend. The horizontal velocity vector field distribution of the monitoring points around the foundation pit is shown below. Figure 3 As shown.

[0085] In an optional embodiment, the step of estimating the velocity gradient tensor based on the three-dimensional velocity field, extracting the principal strain rate representing the degree of deformation based on the symmetric strain rate component of the velocity gradient tensor, and mapping the principal strain rate corresponding to the triangular element to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point includes:

[0086] A Delaunay triangulation was constructed with all GNSS monitoring points as vertices, and triangular units that crossed the cavities inside the pit were removed to obtain the continuous topology around the pit.

[0087] For each triangular element of the Delaunay triangulation, a system of partial differential equations is established based on the planar coordinates of the vertices and the corresponding horizontal velocity components, and the two-dimensional velocity gradient tensor of each element is obtained by solving the equations.

[0088] The symmetric part of the two-dimensional velocity gradient tensor and its transpose matrix is ​​obtained to obtain the two-dimensional strain rate tensor. The eigenvalues ​​of the two-dimensional strain rate tensor are calculated, and the eigenvalue with the largest absolute value is extracted as the principal strain rate of the triangular element.

[0089] The principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are statistically fused to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point.

[0090] Using all GNSS measuring points around the foundation pit, an initial Delaunay triangulation network was constructed based on the Bowyer-Watson point-by-point insertion algorithm. For the non-convex hull case including the central excavation area of ​​the foundation pit, the intersection determination of polygon boundaries was used. When any side length inside the triangulation network is greater than the set threshold of 50 meters, or when any side line intersects geometrically with the boundary of the designed void inside the foundation pit, the triangular unit is removed from the topology, resulting in a continuous topological structure of approximately 120 effective triangles closely surrounding the foundation pit.

[0091] For each selected triangular element, the two-dimensional planar coordinates of its three vertices and the corresponding eastward and northward horizontal velocity components are obtained. Based on the uniform strain assumption of continuum mechanics, a partial derivative of the horizontal velocity with respect to the planar coordinates is constructed, i.e. , , , The system of linear partial differential equations is used. The system of equations is solved using the least squares method or a linear equation solving method to obtain the specific values ​​of the two-dimensional velocity gradient tensor matrix L constituting the triangular element. The calculated two-dimensional velocity gradient tensor matrix L and its transpose are symmetricized to obtain the two-dimensional strain rate tensor matrix. The eigenvalues ​​in two orthogonal directions, namely the first principal strain rate and the second principal strain rate, are obtained by solving the characteristic equation of this matrix. The eigenvalue with the largest absolute value among these two eigenvalues ​​is extracted, for example, a value of 1.5 × 10⁻⁶. -4 / day, this value is set as the principal strain rate scalar of the unit region, as a parameter to quantitatively represent the degree of local differential deformation and shear deformation of the foundation pit.

[0092] For each GNSS monitoring point, all adjacent triangular elements with that monitoring point as the vertex are retrieved; the principal strain rates of each adjacent triangular element are weighted and averaged using the ratio of the area of ​​each adjacent triangular element to the total area of ​​all adjacent triangular elements as the weight, to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point; when a monitoring point is associated with only one valid triangular element, the principal strain rate of that triangular element is directly used as the principal strain rate of the monitoring point.

[0093] In an optional embodiment, the core point discrimination condition of the density clustering algorithm is adjusted by nonlinear mapping from the principal strain rate to the neighborhood search influence factor, and the updated core point discrimination condition is applied to cluster the three-dimensional velocity field to obtain multiple candidate region partitioning schemes, including:

[0094] An exponential decay function with the dimensionless principal strain rate as the independent variable is constructed, and a neighborhood search influence factor between 0 and 1 is calculated.

[0095] The minimum number of contained points preset in the density clustering algorithm is multiplied by the neighborhood search influence factor to obtain the updated core point discrimination threshold. The larger the principal strain rate, the smaller the neighborhood search influence factor, and the lower the corresponding core point discrimination threshold.

[0096] Given a neighborhood radius, the three-dimensional velocity field is subjected to multiple density clustering scans using an updated core point discrimination threshold to obtain multiple candidate region partitioning schemes.

[0097] Min-Max normalization was used to map the principal strain rate sequences of each monitoring point to the [0,1] interval to obtain the dimensionless independent variable x. An exponential decay function of the form f(x) = exp(-k·x) was constructed, where the decay coefficient k was preferably set to 2.0. The independent variable of each monitoring point was substituted into this decay function to calculate the neighborhood search influence factor, so that the influence factor of the high strain region was relatively reduced. For example, when x = 0.8, the calculated influence factor was 0.201.

[0098] Based on the standard DBSCAN density clustering algorithm, a globally preset minimum threshold for the number of included points is set to 5. This threshold is multiplied point-to-point by the neighborhood search influence factor corresponding to each measured point, and then rounded up. This result is then compared with the preset minimum core point threshold, and the larger value is used as the updated core point discrimination threshold. When a measured point is in a region of severe deformation, the core point discrimination threshold is automatically reduced, for example, from 5 to... This makes it easier to identify stations with prominent deformation anomalies but sparse point distribution as density cores, thereby reducing the impact of background noise on the identification of strong deformation features.

[0099] The candidate set of Euclidean distance neighborhood radii in the three-dimensional velocity feature space is set to [0.5, 1.0, 1.5, 2.0] mm / yr. Connectivity expansion clustering is performed using these four neighborhood radii in combination with the aforementioned adjusted dynamic core point discrimination threshold. Each complete full network scan will generate a topology structure containing different cluster combinations, outputting four candidate region partitioning schemes corresponding to these four radii, for subsequent cost function to perform global optimal solution matching.

[0100] S2, for each candidate region partitioning scheme, construct a joint penalty function, traverse all candidate region partitioning schemes and select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, for each region of the optimal region partitioning scheme, according to the mean square error of the observation values ​​of each monitoring point and the model prediction residual, use the iterative reweighted least squares method to solve the local translation parameters and local rotation parameters of each region.

[0101] For each candidate region division scheme, first determine whether each sub-region satisfies the local motion parameter solution conditions; if the number of monitoring points in a certain sub-region is less than the preset minimum number of solution points, or if the design matrix of its local motion model does not meet the preset rank condition, the candidate region division scheme is determined to be an invalid scheme and removed from the candidate scheme set.

[0102] In the candidate region partitioning schemes that meet the solution conditions, monitoring points at the boundaries of adjacent clustering regions are extracted. The predicted velocity difference at the boundary points is calculated based on the initial rigid body motion models of each adjacent region, and the sum of the L2 norms of these predicted velocity differences is used as the residual velocity discontinuity norm. The residual velocity discontinuity norm is then dimensionlessly processed using the baseline velocity norm to obtain a data fitting term. Simultaneously, a model complexity penalty term is constructed based on the current number of region partitions, and a dimensionless reduction weight is generated based on the dimensionless boundary principal strain rate statistics to reduce the model complexity penalty term. The data fitting term is added to the reduced model complexity penalty term to obtain a joint penalty function. All candidate region partitioning schemes are then iterated through, and the scheme with the smallest joint penalty function value is selected as the optimal region partitioning scheme.

[0103] For each sub-region in the optimal region partitioning scheme, local rigid body translation and rotation motion observation equations are established. The reciprocal of the square of the velocity error at each monitoring point is used as the initial observation weight, and the parameters are solved using an iterative reweighted least squares method. In each iteration, a standardized residual is calculated based on the residual between the observed velocity and the model-predicted velocity. A robust weight function is constructed based on the standardized residual, and the weight matrix is ​​updated using the robust weight function. When the changes in local translation parameters and local rotation parameters obtained from two consecutive iterations are less than the corresponding preset convergence thresholds, the local translation parameters and local rotation parameters of that region are output.

[0104] In an optional embodiment, the step of constructing a joint penalty function for each candidate region partitioning scheme, and traversing all candidate region partitioning schemes to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, includes:

[0105] The validity of the sub-regions in each candidate region partitioning scheme is determined, and candidate region partitioning schemes that do not meet the preset minimum number of solution points or preset rank conditions are eliminated.

[0106] The residual velocity discontinuity norm of all monitoring points on the boundary of adjacent regions in each candidate region division scheme is calculated as a data fitting term;

[0107] A model complexity penalty term is constructed based on the number of regions. The boundary principal strain rate statistics are calculated according to the principal strain rate of the monitoring points corresponding to the boundary monitoring points. The boundary principal strain rate statistics are then dimensionless to generate dimensionless reduction weights. These dimensionless reduction weights are then used to reduce the model complexity penalty term.

[0108] The data fitting term after being dimensionless by the benchmark velocity norm is added to the reduced dimensionless model complexity penalty term to form the joint penalty function. All candidate region partitioning schemes are traversed, and the scheme with the smallest value of the joint penalty function is selected as the optimal region partitioning scheme. The number of regions at this time is the optimal number of local motion models.

[0109] For any clustering partitioning scheme, first check whether each independent region meets the local motion parameter solution conditions. If the number of monitoring points in any sub-region is less than the preset minimum number of solution points, or the design matrix of the corresponding local motion model does not meet the preset rank condition, then the candidate region partitioning scheme is determined to be an invalid scheme and removed from the candidate scheme set. For the remaining valid candidate region partitioning schemes, calculate the residual velocity discontinuity norm and model complexity penalty term.

[0110] For candidate region partitioning schemes that meet the solution conditions, firstly, for each sub-region, an initial local rigid body motion model is fitted using the three-dimensional coordinates and three-dimensional velocity observations of the monitoring points within that sub-region; then, boundary monitoring point pairs belonging to adjacent sub-regions and connected by Delaunay edges are retrieved, and the spatial midpoint of the boundary monitoring point pair is calculated as the boundary representative point; using the initial local rigid body motion models of the two adjacent sub-regions respectively, the three-dimensional velocity vectors are predicted at the boundary representative point, and the Euclidean norm 2 of the difference between the two three-dimensional predicted velocity vectors is calculated as the velocity discontinuity of the boundary monitoring point pair. The norm 2 of all cross-region boundary point pairs in the entire network is scalarly summed, for example, accumulating to 15.2 mm / yr. This summation value is used as a data fitting term to evaluate the degree of velocity jump at the model splicing point.

[0111] Extract the total number of independent regions divided by this scheme. Input the total number of points Construct the initial complexity penalty formula The principal strain rates corresponding to all boundary monitoring points are extracted and their average values ​​are calculated. These average principal strain rates are then normalized using a preset reference strain rate to obtain dimensionless boundary principal strain rate statistics. Design the attenuation weight formula , where β is a dimensionless adjustment coefficient. The calculated dimensionless weight factor w is multiplied by the initial penalty formula to reduce the model complexity penalty term, so that the system tends to retain necessary partitions in the strong deformation zone, reducing the risk of underpartition caused by excessive complexity penalty.

[0112] The root mean square of the actual observed 3D velocities across the entire network is extracted as the baseline velocity, such as 3.5 mm / yr. The data fitting term is divided by this baseline velocity to achieve dimensionless processing. The dimensionless fitting term is then added to the previously reduced model complexity penalty term to establish a joint penalty function formula. The four candidate partitioning schemes generated in the previous steps are evaluated, yielding penalty function values ​​of [4.2, 3.5, 2.1, 5.0]. Through numerical comparison, the partitioning scheme with a penalty value of 2.1 is selected as the optimal region partitioning scheme. Furthermore, the number of regions included in this scheme, for example, dividing the area into three independent motion blocks, is considered the optimal number of local motion models for the current deformation analysis.

[0113] In an optional embodiment, the step of calculating the local translation and local rotation parameters of each region using iterative reweighted least squares based on the mean square error of the observed values ​​at each monitoring point and the model prediction residual includes:

[0114] An initial diagonal weight matrix is ​​constructed using the reciprocal of the squared error of the observation values ​​at each monitoring point, and the initial local translation parameters and initial local rotation parameters are initially calculated.

[0115] The model-predicted velocity is calculated based on the initial local translation parameters and the initial local rotation parameters, and the vector residual between the model-predicted velocity and the actual observed velocity is obtained.

[0116] The standardized residuals are calculated based on the vector residuals and the post-correction unit weight mean square error of the overall adjustment, and a robust weight function is constructed based on the standardized residuals to obtain a new weight matrix.

[0117] The new weight matrix is ​​used to replace the original weight matrix for iterative solution until the change in local translation parameters between two consecutive solutions is less than the preset translation convergence threshold, and the change in local rotation parameters between two consecutive solutions is less than the preset rotation convergence threshold. The local translation parameters and local rotation parameters are then output.

[0118] For a specific block in the optimal region partitioning scheme, extract all monitoring points contained within it, for example, 15 monitoring points; read the estimated mean square error (MSE) of the northward, eastward, and vertical velocity components of each monitoring point from the three-dimensional velocity field calculation results, for example, 0.02 mm / yr; use the square of the estimated mean square error as the velocity observation variance, and use the reciprocal of the velocity observation variance as the weight to fill and construct an initial diagonal weight matrix. A Jacobian matrix A is constructed that maps the observed velocity to six unknown motion parameters. The classical weighted least squares formula is then applied. Solving for the included Three translations and The initial 6-parameter model vector with three rotations.

[0119] Substituting the calculated initial 6-parameter vector into the local rigid body motion model, the three-dimensional predicted velocity vectors of each monitoring point within the region are calculated. The actual observed velocity is then subtracted from the predicted velocity, and the Euclidean distance is taken to calculate the independent velocity vector residual for each point. Using the IGG III robust estimation criterion, the apologetic unit weighted mean square error of the overall adjustment (e.g., 0.5 mm / yr) is input to calculate the standardized residual for each measuring point. A piecewise robust weight function is constructed based on the standardized residual. When the absolute value of the standardized residual is not greater than the reduction threshold... When the weights remain unchanged, the absolute value of the standardized residuals is greater than the reduction threshold. And not greater than the elimination threshold At that time, according to The weights decay exponentially. Indicates the first The absolute value of the standardized residuals at each measurement point; when the absolute value of the standardized residuals exceeds the elimination threshold. When this is done, the weight is set to a preset minimum value ε. Preferably, =1.5, =3.0, ε=10 -6 A new diagonal weight matrix is ​​constructed. .

[0120] Using the new diagonal weight matrix Replace the equation The solution process is repeated, iterating through the steps of predictive velocity calculation, residual calculation, and weight matrix update. After each solution, the changes in the translation parameter vector and rotation parameter vector obtained in the current and previous iterations are calculated, and it is determined whether they are less than the corresponding preset translation convergence threshold and preset rotation convergence threshold, respectively. If, after multiple iterations, such as four iterations, the changes in the translation and rotation parameters both satisfy the corresponding convergence conditions, the iteration ends, and the state vector at the point of convergence is output as the local translation and local rotation parameters for that independent region.

[0121] S3. Perform topology optimization iteration operation based on the local translation parameters and local rotation parameters of each region. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency condition, merge adjacent regions and reassign the boundary monitoring points with smaller fitting residuals to adjacent regions. Alternately execute region merging and boundary monitoring point reassignment until the region division is stable, and output the foundation pit deformation field model.

[0122] A regional adjacency graph is constructed based on the regional spatial adjacency relationships in the optimal regional partitioning scheme. Nodes represent sub-regions, and edges indicate the adjacency boundaries between two sub-regions. For each edge in the regional adjacency graph, the local translation parameter difference norm and local rotation parameter difference norm between the two adjacent regions connected by it are calculated. If both are lower than preset translation consistency thresholds and preset rotation consistency thresholds, the two adjacent regions are determined to meet the motion consistency condition and are merged into a single region. Then, the local translation parameters and local rotation parameters are recalculated based on the merged set of monitoring points.

[0123] After merging the regions, monitoring points at the boundaries of each region are extracted. Predicted velocities are calculated using the local rigid body motion models of the current region and its adjacent regions for each boundary monitoring point. The difference between the predicted velocity and the actual observed velocity is used as the fitting residual. If the fitting residual of a boundary monitoring point relative to its adjacent regions is less than its fitting residual relative to its current region, and the reduction ratio meets a preset residual reduction condition, then the boundary monitoring point's label is updated, and it is reassigned to its corresponding adjacent region.

[0124] The system alternately performs region merging and boundary monitoring point reassignment operations, and determines whether the number of regions and the monitoring point affiliation labels have changed after each iteration. If there are no adjacent regions that meet the motion consistency condition during continuous iterations, and no boundary monitoring points have changed their affiliation, the region division is determined to have reached a stable state.

[0125] The converged boundary information of each region, along with the corresponding local translation and rotation parameters, are combined and encapsulated to generate a spatially discrete kinematic parameter field data structure. This structure serves as the foundation pit deformation field model for visualization output and subsequent disaster early warning.

[0126] In an optional embodiment, the step of performing topology optimization iteration based on the local translation parameters and local rotation parameters of each region, merging adjacent regions when the differences in translation and rotation parameters between adjacent regions meet a preset motion consistency condition, and reassigning boundary monitoring points with smaller fitting residuals to adjacent regions, alternately performing region merging and boundary monitoring point reassignment until the region division is stable, includes:

[0127] Calculate the difference in local translation vectors and the difference in local rotation vectors between two adjacent regions. Based on the average enclosing radius of the two adjacent regions, convert the difference in local rotation vectors into an equivalent linear velocity difference. Calculate the composite difference norm based on the difference in local translation vectors and the equivalent linear velocity difference.

[0128] If the composite difference norm is lower than the preset composite consistency threshold, then the two adjacent regions are merged, and the local translation parameters and local rotation parameters of the merged region are recalculated.

[0129] Calculate the fitting residual of the boundary monitoring point to the local motion model of the original region to which the monitoring point belongs, and the fitting residual of the local motion model of all adjacent regions. If the fitting residual of a certain adjacent region is less than the fitting residual of the original region, and the difference between the two satisfies the preset residual reduction condition, then the monitoring point is reclassified to the adjacent region.

[0130] The process alternately executes the region merging and boundary monitoring point reassignment until no adjacent region meets the preset motion consistency condition and no monitoring point undergoes a change in assignment. Then, the iteration terminates and the foundation pit deformation field model is output.

[0131] For any two adjacent regions with spatially connected boundaries and The differences in local translation vectors and local rotation vectors of the two regions are calculated separately. Using the average enclosing radius of the two regions as a scale factor, the difference in local rotation vectors is converted into an equivalent linear velocity difference. Then, the sum of the squares of the local translation vector difference and the equivalent linear velocity difference is taken, and the square root is obtained to obtain a dimensionlessly consistent composite difference norm. A consistency motion judgment threshold of 2.0 mm / yr is set. If the calculated difference norm is 1.2 mm / yr, which is lower than this threshold, then... and The GNSS measurement point list is merged, and the robust least squares module is called to refit a 6-parameter local rigid body motion model for the merged single region.

[0132] Extract all boundary monitoring points. Substitute the spatial coordinates of each boundary monitoring point into the local rigid body motion model of its original region and adjacent candidate regions, respectively, calculate the corresponding predicted velocity, and compare the predicted velocity with the observed three-dimensional velocity of the boundary monitoring point to obtain the model fitting residuals of the original region and the adjacent region. If the model fitting residual of the adjacent region is less than that of the original region, and the reduction ratio of the two meets the preset residual reduction condition, then the measuring point is reassigned to the adjacent deformation region with better performance. The region merging operation based on motion consistency and the boundary point reassignment operation based on residual reduction condition are executed alternately. If there is region fusion or cross-boundary migration of measuring points, the parameters of the entire network are updated and the next iteration begins. The above scanning and detection steps are executed alternately until no region with a difference norm lower than 2.0 mm / yr is detected and no measuring point transfer that meets the residual reduction condition occurs. Then, the iteration stops, and the partition status and corresponding parameters at this time are output as the foundation pit deformation field model.

[0133] The experiment selected 120 continuous monitoring points deployed around a large deep foundation pit project as the test dataset, with a monitoring time span of 180 days, acquiring complete three-dimensional coordinate time series records. The experimental hardware was equipped with a high-performance computing node with a processor clock speed of 3.2 GHz. To verify the effectiveness of the proposed method, three sets of control experiments were designed. Scheme 1 is the traditional basic model, which only uses the standard density clustering algorithm to fix the core point discrimination threshold at 5 and uses ordinary least squares method for parameter calculation. Scheme 2, based on Scheme 1, adjusts the clustering core point threshold by using the nonlinear mapping from principal strain rate to neighborhood search influence factors. Scheme 3 is the optimal model fully proposed in this invention, which uses iterative reweighted least squares method and robust weight function for parameter refinement calculation based on Scheme 2.

[0134] Regarding the accuracy of region segmentation, Scheme 1, due to its fixed threshold, resulted in numerous missed detections in areas with severe deformation but sparse data points, achieving an accuracy of only 76.5%. Scheme 2, by adjusting the threshold for high-strain zones using an exponential decay function, improved the ability to identify strongly deformed zones, increasing the region segmentation accuracy to 89.2%. In terms of the root mean square error (RMSE) of velocity field prediction, Scheme 1 was significantly affected by gross errors and outliers, reaching an error of 2.8 mm / yr. Scheme 2 reduced the error to 2.4 mm / yr. Scheme 3, by constructing a new diagonal weight matrix using the post-hoc unit weight mean square error, weakened the impact of gross errors, reducing the RMS error of velocity field prediction to 1.1 mm / yr, and meeting the accuracy requirement of less than the set convergence condition within four consecutive iterations.

[0135] The improvement of Scheme 2 compared to Scheme 1 confirms the effectiveness of the core point discrimination condition. Constructing dimensionless independent variables based on principal strain rates and generating influencing factors using attenuation coefficients makes monitoring stations with prominent deformation anomalies more easily identified as density cores, reducing the impact of background noise on the identification of strong deformation features. The improvement of Scheme 3 compared to Scheme 2 is attributed to the application of the iterative reweighted least squares criterion. By alternately performing residual calculations and weight matrix updates, the weights of outliers with excessively large standardized residuals are attenuated or even set to minimum values. This avoids the local translation and rotation parameter calculation process being driven by local distortion points, and the output model can better represent the observed overall surface movement trend.

[0136] In the second embodiment, the present invention also proposes a foundation pit deformation monitoring system, comprising the following modules:

[0137] The processing module is used to acquire the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit, calculate the three-dimensional velocity field of each monitoring point, estimate the velocity gradient tensor based on the three-dimensional velocity field, extract the principal strain rate representing the degree of deformation according to the symmetric strain rate component of the velocity gradient tensor, map the principal strain rate corresponding to the triangular unit to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point, adjust the core point discrimination condition of the density clustering algorithm through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor, apply the updated core point discrimination condition to cluster the three-dimensional velocity field, and obtain multiple candidate region division schemes.

[0138] The module is used to construct a joint penalty function for each candidate region partitioning scheme, traverse all candidate region partitioning schemes and select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, and for each region of the optimal region partitioning scheme, calculate the local translation parameters and local rotation parameters of each region based on the mean square error of the observation values ​​of each monitoring point and the model prediction residual using the iterative reweighted least squares method.

[0139] The output module is used to perform topology optimization iteration based on the local translation parameters and local rotation parameters of each region. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency conditions, the adjacent regions are merged. The boundary monitoring points with smaller fitting residuals to the adjacent regions are reassigned. The region merging and boundary monitoring point reassignment are alternately executed until the region division is stable, and the foundation pit deformation field model is output.

[0140] It should be understood that the various forms of processes shown above can be used to rearrange, add, or delete steps. For example, the steps described in this disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this disclosure can be achieved, and this is not limited herein.

[0141] The specific embodiments described above do not constitute a limitation on the scope of protection of this disclosure. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the principles of this disclosure should be included within the scope of protection of this disclosure.

Claims

1. A method for monitoring the deformation of a foundation pit, characterized in that, include: The three-dimensional coordinate time series of GNSS monitoring points around the foundation pit is obtained, and the three-dimensional velocity field of each monitoring point is calculated. The velocity gradient tensor is estimated based on the three-dimensional velocity field, and the principal strain rate representing the degree of deformation is extracted according to the symmetric strain rate component of the velocity gradient tensor. The principal strain rate corresponding to the triangular unit is mapped to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point. Through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor, the core point discrimination condition of the density clustering algorithm is adjusted. The updated core point discrimination condition is applied to the three-dimensional velocity field for clustering processing to obtain multiple candidate region division schemes. For each candidate region partitioning scheme, a joint penalty function is constructed. All candidate region partitioning schemes are traversed to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme. For each region of the optimal region partitioning scheme, the local translation parameters and local rotation parameters of each region are calculated using the iterative reweighted least squares method based on the mean square error of the observation values ​​of each monitoring point and the model prediction residual. Based on the local translation parameters and local rotation parameters of each region, a topology optimization iterative operation is performed. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency condition, the adjacent regions are merged. The boundary monitoring points with smaller fitting residuals for adjacent regions are reassigned. The region merging and boundary monitoring point reassignment are alternately executed until the region division is stable, and the foundation pit deformation field model is output.

2. The method according to claim 1, characterized in that, The process of obtaining the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit and calculating the three-dimensional velocity field of each monitoring point includes: The three-dimensional coordinate time series is denoised and gross error removed to obtain a smoothed three-dimensional coordinate time series. The smoothed three-dimensional coordinate time series was linearly fitted using the sliding window method, and the deformation rate within each time window was calculated. The deformation rate within each time window is time-weighted and averaged to obtain the three-dimensional velocity of each monitoring point in the north, east and vertical directions, and combined to obtain the three-dimensional velocity field.

3. The method according to claim 1, characterized in that, The process of estimating the velocity gradient tensor based on the three-dimensional velocity field, extracting the principal strain rates representing the degree of deformation based on the symmetric strain rate components of the velocity gradient tensor, and mapping the principal strain rates corresponding to the triangular elements to the principal strain rates of the monitoring points corresponding to the GNSS monitoring points includes: A Delaunay triangulation was constructed with all GNSS monitoring points as vertices, and triangular units that crossed the cavities inside the pit were removed to obtain the continuous topology around the pit. For each triangular element of the Delaunay triangulation, a system of partial differential equations is established based on the planar coordinates of the vertices and the corresponding horizontal velocity components, and the two-dimensional velocity gradient tensor of each element is obtained by solving the equations. The symmetric part of the two-dimensional velocity gradient tensor and its transpose matrix is ​​obtained to obtain the two-dimensional strain rate tensor. The eigenvalues ​​of the two-dimensional strain rate tensor are calculated, and the eigenvalue with the largest absolute value is extracted as the principal strain rate of the triangular element. The principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are statistically fused to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point.

4. The method according to claim 2 or 3, characterized in that, The nonlinear mapping from the principal strain rate to the neighborhood search influence factor is used to adjust the core point discrimination condition of the density clustering algorithm. The updated core point discrimination condition is then applied to cluster the three-dimensional velocity field, resulting in multiple candidate region partitioning schemes, including: An exponential decay function with the dimensionless principal strain rate as the independent variable is constructed, and a neighborhood search influence factor between 0 and 1 is calculated. The minimum number of contained points preset in the density clustering algorithm is multiplied by the neighborhood search influence factor to obtain the updated core point discrimination threshold. The larger the principal strain rate, the smaller the neighborhood search influence factor, and the lower the corresponding core point discrimination threshold. Given a neighborhood radius, the three-dimensional velocity field is subjected to multiple density clustering scans using an updated core point discrimination threshold to obtain multiple candidate region partitioning schemes.

5. The method according to claim 1, characterized in that, The process of constructing a joint penalty function for each candidate region partitioning scheme, and then iterating through all candidate region partitioning schemes to select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, includes: The validity of the sub-regions in each candidate region partitioning scheme is determined, and candidate region partitioning schemes that do not meet the preset minimum number of solution points or preset rank conditions are eliminated. The residual velocity discontinuity norm of all monitoring points on the boundary of adjacent regions in each candidate region division scheme is calculated as a data fitting term; A model complexity penalty term is constructed based on the number of regions. The boundary principal strain rate statistics are calculated according to the principal strain rate of the monitoring points corresponding to the boundary monitoring points. The boundary principal strain rate statistics are then dimensionless to generate dimensionless reduction weights. These dimensionless reduction weights are then used to reduce the model complexity penalty term. The data fitting term after being dimensionless by the benchmark velocity norm is added to the reduced dimensionless model complexity penalty term to form the joint penalty function. All candidate region partitioning schemes are traversed, and the scheme with the smallest value of the joint penalty function is selected as the optimal region partitioning scheme. The number of regions at this time is the optimal number of local motion models.

6. The method according to claim 1, characterized in that, The method involves calculating the local translation and rotation parameters of each region using iterative reweighted least squares based on the mean square error of the observed values ​​at each monitoring point and the model prediction residual. This includes: An initial diagonal weight matrix is ​​constructed using the reciprocal of the squared error of the observation values ​​at each monitoring point, and the initial local translation parameters and initial local rotation parameters are initially calculated. The model-predicted velocity is calculated based on the initial local translation parameters and the initial local rotation parameters, and the vector residual between the model-predicted velocity and the actual observed velocity is obtained. The standardized residuals are calculated based on the vector residuals and the post-correction unit weight mean square error of the overall adjustment, and a robust weight function is constructed based on the standardized residuals to obtain a new weight matrix. The new weight matrix is ​​used to replace the original weight matrix for iterative solution until the change in local translation parameters between two consecutive solutions is less than the preset translation convergence threshold, and the change in local rotation parameters between two consecutive solutions is less than the preset rotation convergence threshold. The local translation parameters and local rotation parameters are then output.

7. The method according to claim 1, characterized in that, The topology optimization iterative operation based on the local translation parameters and local rotation parameters of each region, merging adjacent regions when the differences in translation and rotation parameters between adjacent regions meet a preset motion consistency condition, and reassigning the boundary monitoring points with smaller fitting residuals for adjacent regions, alternately performing region merging and boundary monitoring point reassignment until the region division is stable, includes: Calculate the difference in local translation vectors and the difference in local rotation vectors between two adjacent regions. Based on the average enclosing radius of the two adjacent regions, convert the difference in local rotation vectors into an equivalent linear velocity difference. Calculate the composite difference norm based on the difference in local translation vectors and the equivalent linear velocity difference. If the composite difference norm is lower than the preset composite consistency threshold, then the two adjacent regions are merged, and the local translation parameters and local rotation parameters of the merged region are recalculated. Calculate the fitting residual of the boundary monitoring point to the local motion model of the original region to which the monitoring point belongs, and the fitting residual of the local motion model of all adjacent regions. If the fitting residual of a certain adjacent region is less than the fitting residual of the original region, and the difference between the two satisfies the preset residual reduction condition, then the monitoring point is reclassified to the adjacent region. The process alternately executes the region merging and boundary monitoring point reassignment until no adjacent region meets the preset motion consistency condition and no monitoring point undergoes a change in assignment. Then, the iteration terminates and the foundation pit deformation field model is output.

8. A foundation pit deformation monitoring system, characterized in that, include: The processing module is used to acquire the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit, calculate the three-dimensional velocity field of each monitoring point, estimate the velocity gradient tensor based on the three-dimensional velocity field, extract the principal strain rate representing the degree of deformation according to the symmetric strain rate component of the velocity gradient tensor, map the principal strain rate corresponding to the triangular unit to the principal strain rate of the monitoring point corresponding to the GNSS monitoring point, adjust the core point discrimination condition of the density clustering algorithm through the nonlinear mapping from the principal strain rate to the neighborhood search influence factor, apply the updated core point discrimination condition to cluster the three-dimensional velocity field, and obtain multiple candidate region division schemes. The module is used to construct a joint penalty function for each candidate region partitioning scheme, traverse all candidate region partitioning schemes and select the scheme that minimizes the value of the joint penalty function as the optimal region partitioning scheme, and for each region of the optimal region partitioning scheme, calculate the local translation parameters and local rotation parameters of each region based on the mean square error of the observation values ​​of each monitoring point and the model prediction residual using the iterative reweighted least squares method. The output module is used to perform topology optimization iteration based on the local translation parameters and local rotation parameters of each region. When the differences in translation parameters and rotation parameters between adjacent regions meet the preset motion consistency conditions, the adjacent regions are merged. The boundary monitoring points with smaller fitting residuals to the adjacent regions are reassigned. The region merging and boundary monitoring point reassignment are alternately executed until the region division is stable, and the foundation pit deformation field model is output.

9. The system according to claim 8, characterized in that, The process of obtaining the three-dimensional coordinate time series of GNSS monitoring points around the foundation pit and calculating the three-dimensional velocity field of each monitoring point includes: The three-dimensional coordinate time series is denoised and gross error removed to obtain a smoothed three-dimensional coordinate time series. The smoothed three-dimensional coordinate time series was linearly fitted using the sliding window method, and the deformation rate within each time window was calculated. The deformation rate within each time window is time-weighted and averaged to obtain the three-dimensional velocity of each monitoring point in the north, east and vertical directions, and combined to obtain the three-dimensional velocity field.

10. The system according to claim 8, characterized in that, The process of estimating the velocity gradient tensor based on the three-dimensional velocity field, extracting the principal strain rates representing the degree of deformation based on the symmetric strain rate components of the velocity gradient tensor, and mapping the principal strain rates corresponding to the triangular elements to the principal strain rates of the monitoring points corresponding to the GNSS monitoring points includes: A Delaunay triangulation was constructed with all GNSS monitoring points as vertices, and triangular units that crossed the cavities inside the pit were removed to obtain the continuous topology around the pit. For each triangular element of the Delaunay triangulation, a system of partial differential equations is established based on the planar coordinates of the vertices and the corresponding horizontal velocity components, and the two-dimensional velocity gradient tensor of each element is obtained by solving the equations. The symmetric part of the two-dimensional velocity gradient tensor and its transpose matrix is ​​obtained to obtain the two-dimensional strain rate tensor. The eigenvalues ​​of the two-dimensional strain rate tensor are calculated, and the eigenvalue with the largest absolute value is extracted as the principal strain rate of the triangular element. The principal strain rates of multiple triangular elements adjacent to the same GNSS monitoring point are statistically fused to obtain the principal strain rate of the monitoring point corresponding to that GNSS monitoring point.

Images