A foundation pit deformation space-time evolution prediction and support parameter reverse optimization system and method

Abstract

The application relates to the technical field of geotechnical engineering and foundation pit supporting design, and particularly discloses a system and method for predicting spatiotemporal evolution of foundation pit deformation and reversely optimizing supporting parameters, which comprises the following steps: collecting sparse borehole geotechnical parameters and in-situ test data to construct a prior random field; adopting Bayesian compressed sensing to invert parameter space full-grid mean values and covariance matrices under posterior probability distribution; generating multiple groups of random fields to realize and input finite element simulation, forming a training sample set of parameter distribution and deformation field; constructing a deep Gaussian process nonlinear mapping and quantifying prediction variance; and reversely searching supporting parameter combinations by taking deformation control threshold and allowable overrun probability as constraints. The application can quantitatively influence the variability of geological parameter space, and obtain a robust supporting scheme meeting overrun probability requirements.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering and foundation pit support design technology, specifically to a system and method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters. Background Technology

[0002] During the excavation process of foundation pit engineering, soil deformation exhibits significant spatiotemporal evolution characteristics. Reasonable deformation prediction and subsequent optimization of support parameters are crucial to ensuring the safety of the foundation pit and its surrounding environment. Existing technologies commonly used methods for predicting foundation pit deformation include finite element numerical simulation, empirical formula methods, and time-series prediction models based on monitoring data.

[0003] The existing technology has the following shortcomings: Under complex geological conditions with significant spatial variability in soil and rock parameters, existing deterministic inverse optimization methods cannot quantify the impact of geological uncertainties on deformation response, resulting in an uncontrollable probability of deformation exceeding limits when the optimized support parameter combination is exposed to local anomalous geological bodies during actual construction. Summary of the Invention

[0004] The purpose of this invention is to provide a system and method for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters, so as to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A method for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters includes the following steps: S1. Collect sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area to construct a prior random field containing parameter mean, variance and autocorrelation distance; S2 uses Bayesian compressed sensing to perform sparse observation inversion on the prior random field and outputs the mean and covariance matrix of the parameter space full grid under the posterior probability distribution. S3, based on the posterior probability distribution, generates multiple sets of random fields. Each set of random fields is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. S4. Construct a nonlinear mapping function from the parameter space distribution to the deformable field using the training sample set, and quantify the prediction variance of the mapping output. S5 uses the deformation control threshold and the allowable over-limit probability as constraints to inversely search for the support parameter combination, so that the probability of satisfying the constraint conditions under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

[0006] As a further aspect of the present invention: S1 specifically includes: The sparse borehole geotechnical parameters are used as hard data, and the continuous in-situ test data are used as soft data. The initial distribution of the parameter mean is generated by fusion. The anisotropic variability function is calculated based on the parameter mean distribution, and the spatial autocorrelation distance is then derived. The parameter variance is locally corrected using cross-validation residuals, and a prior random field is output.

[0007] As a further aspect of the present invention: S2 specifically includes: Discretize the autocorrelation matrix of the prior random field into a sparse basis matrix and establish a linear mapping relationship between the observed data and the grid parameters; We introduce Laplace prior constraints and use variational Bayesian inference to alternately update hyperparameters and posterior covariance matrices. Stop when the change in the posterior mean between two consecutive iterations is less than a set threshold, and output the full grid mean and covariance matrix of the parameter space under the posterior probability distribution.

[0008] As a further aspect of the present invention: the introduction of Laplace prior constraints and the use of variational Bayesian inference to alternately update hyperparameters and posterior covariance matrices specifically include: A sparse linear observation equation is constructed based on the autocorrelation matrix of the prior random field and the observation data, and a Laplace prior distribution is assigned to the full grid mean of the parameter space. The Laplace prior distribution is decomposed into a Gaussian scale mixture, and an auxiliary scale parameter is introduced and its gamma prior is set. The variational Bayesian expectation-maximization step is adopted to sequentially update the posterior expectation, posterior covariance matrix and hyperparameter of the auxiliary scaling parameter until the difference between the posterior means of two adjacent iterations is less than the preset tolerance, and then the updated posterior mean and covariance matrix are output.

[0009] As a further aspect of the present invention: S3 specifically includes: The covariance matrix under the posterior probability distribution is decomposed by Cholliski decomposition, and combined with the full grid mean, a predetermined number of random fields are generated by sequential Gaussian simulation. For each set of random fields, the corresponding excavation sequence and support conditions are matched, and the finite element solver is called sequentially to calculate the spatial deformation field and time history curve of each construction stage. The grid node values ​​of the parameter space distribution are stored one pair at a time with the corresponding deformation field node values ​​to form a training sample set.

[0010] As a further aspect of the present invention: the generation of a predetermined number of random fields specifically includes: The posterior covariance matrix is ​​subjected to incomplete Choleski decomposition to obtain a low-rank approximation lower triangular decomposition matrix. Conditional simulations are performed sequentially across all grid nodes. After each node's parameter value is generated, the covariance conditions of adjacent unsampled nodes are immediately adjusted based on the updated values ​​of the corresponding nodes. Repeat the sequential process until all grid nodes are assigned values, and output a set of random field implementations that follow a posterior probability distribution.

[0011] As a further aspect of the present invention: S4 specifically includes: The parameter spatial distribution in the training sample set is used as input and the deformation field is used as output to construct a hierarchical structure of the deep Gaussian process. Each layer adopts the induced point sparsity approximation to reduce the computational complexity. The length scale, signal variance, and induced point location of the kernel function of each layer are jointly learned by variational Bayesian inference, while the variational lower bound of the marginal likelihood is optimized. Based on the learned deep Gaussian process, the posterior prediction distribution is calculated for any parameter distribution in the input space, and the predicted mean and prediction covariance matrix of the output deformation field are obtained. The diagonal elements of the covariance matrix are extracted as the prediction variance.

[0012] As a further aspect of the present invention: the method of jointly learning the length scale, signal variance, and induced point location of each layer kernel function through variational Bayesian inference, while simultaneously optimizing the lower bound of the marginal likelihood variation, specifically includes: The marginal likelihood variational lower bound is decomposed into the sum of local lower bounds on each training subset, and the gradient of the lower bound with respect to the variational parameters is estimated by stochastic natural gradient descent. With the kernel function's length scale and signal variance fixed, the Riemann gradient ascent update is performed on the induced point position based on the current variational posterior covariance matrix. Alternately update the kernel function hyperparameters and induced point positions until the relative increment of the variational lower bound between two adjacent iterations is lower than a preset threshold, and output the kernel function parameters and induced point positions of each layer after joint learning.

[0013] As a further aspect of the present invention: S5 specifically includes: Based on the deformation mean and prediction variance output by the nonlinear mapping function, an expected improved acquisition function is constructed, which includes deformation control threshold constraints and over-limit penalties. An adaptive evolution strategy based on the covariance matrix is ​​used to iteratively search within the support parameter space. In each iteration, the acquisition function is called to evaluate the deformation over-limit probability corresponding to the candidate parameter combination, and the penalty weight is dynamically adjusted based on the difference between the over-limit probability and the preset level. The search terminates when the relative change of the optimal value of the acquisition function is less than the set tolerance in multiple consecutive iterations, and outputs the combination of support parameters that meets the over-limit probability requirement and minimizes the support cost.

[0014] A system for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters includes: The data acquisition and prior modeling module collects sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area, and constructs a prior random field including parameter mean, variance and autocorrelation distance. The Bayesian compressed sensing inversion module uses Bayesian compressed sensing to perform sparse observation inversion on a prior random field and outputs the mean and covariance matrix of the parameter space under the posterior probability distribution. The random field generation and numerical simulation module generates multiple sets of random field implementations based on the posterior probability distribution. Each set of random field implementations is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. The nonlinear mapping learning module uses the training sample set to construct a nonlinear mapping function from the parameter space distribution to the deformation field, and quantifies the prediction variance of the mapping output. The robust inverse optimization module uses deformation control threshold and allowable over-limit probability as constraints to inversely search for support parameter combinations, so that the probability of satisfying the constraints under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

[0015] The beneficial effects of this invention are: (1) This invention obtains the full-grid posterior probability distribution of the parameter space by constructing a prior random field and using Bayesian compressed sensing for sparse observation inversion. Based on this, multiple random fields are generated to realize and establish a nonlinear mapping function between the parameter space distribution and the deformation field, while quantifying the prediction variance. This method can effectively characterize the spatial variability of soil parameters and transfer geological uncertainty to deformation prediction, so that the support parameter combination obtained by the final inverse optimization can still meet the deformation control threshold with a higher probability than the preset probability under abnormal geological conditions such as local weak interlayers, thereby improving the adaptability of the support scheme to complex strata and the engineering safety.

[0016] (2) This invention utilizes a deep Gaussian process to construct a proxy mapping from the parameter space distribution to the deformation field, and uses variational Bayesian inference to jointly learn the kernel function parameters and induced point locations, significantly reducing the number of calls to finite element numerical simulations. Simultaneously, it employs an adaptive evolution strategy based on the covariance matrix and opportunistic constraint programming for inverse search, automatically finding the optimal solution while ensuring the deformation exceedance probability is controllable. This method avoids the inefficiency of repeatedly calling numerical simulations in traditional trial-and-error or deterministic optimization methods, shortens the support parameter design cycle, and provides confidence intervals for deformation prediction, thus providing technical support for risk quantification and dynamic design in foundation pit engineering. Attached Figure Description

[0017] The invention will now be further described with reference to the accompanying drawings.

[0018] Figure 1This is a flowchart of the method of the present invention; Figure 2 This is a system block diagram of the present invention. Detailed Implementation

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

[0020] Please see Figure 1 As shown, this invention is a method for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters, comprising the following steps: S1. Collect sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area to construct a prior random field containing parameter mean, variance and autocorrelation distance; S2 uses Bayesian compressed sensing to perform sparse observation inversion on the prior random field and outputs the mean and covariance matrix of the parameter space full grid under the posterior probability distribution. S3, based on the posterior probability distribution, generates multiple sets of random fields. Each set of random fields is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. S4. Construct a nonlinear mapping function from the parameter space distribution to the deformable field using the training sample set, and quantify the prediction variance of the mapping output. S5 uses the deformation control threshold and the allowable over-limit probability as constraints to inversely search for the support parameter combination, so that the probability of satisfying the constraint conditions under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

[0021] In S1, sparse borehole geotechnical parameters and continuous in-situ test data are collected in the foundation pit area to construct a prior random field containing parameter mean, variance, and autocorrelation distance, specifically including: First, geotechnical parameters were collected from sparse boreholes in the excavation area. Boreholes were arranged within the excavation area according to exploration specifications, with a spacing of 15 to 30 meters between boreholes. The borehole depth should exceed the excavation depth of the excavation by 1.5 to 2 times. At different depths in each borehole (one sample every 1 to 2 meters), undisturbed soil samples were collected on-site using a thin-walled sampler, sealed, and sent to the indoor geotechnical laboratory. The elastic modulus and cohesion of the soil were obtained through triaxial compression tests, the internal friction angle through direct shear tests, and the compression modulus and void ratio through consolidation tests. The geotechnical parameter values ​​at different depths at each borehole location were obtained. Simultaneously, continuous in-situ testing was conducted in the foundation pit area, specifically using static cone penetration tests: a double-bridge probe (10-15 cm²) connected to a hydraulic static cone penetrometer was used, penetrating downwards at a constant penetration rate of 1-2 cm / s. The probe tip resistance, sidewall friction, and friction ratio were continuously recorded. A set of data was automatically collected every 0.1 meters of penetration, obtaining a continuous parameter profile from the surface to the bottom of the borehole. Standard penetration tests used a 63.5 kg hammer with a drop height of 76 cm, recording the number of blows every 30 cm of penetration, with additional testing at critical depths (such as the soil-rock interface). The aforementioned borehole geotechnical parameters are considered hard data, while the continuous in-situ test data are considered soft data.

[0022] Secondly, the hard and soft data are fused to generate an initial distribution of the parameter mean. The Bayesian maximum entropy method is used for this fusion. Specifically, the hard data is treated as precise observations, while the soft data is transformed into uncertain estimates using a linear regression method. The transformation coefficient is determined by the ratio of the hard data at the same borehole location to the soft data at the corresponding depth, and this ratio is the average of all borehole locations. Then, the entire pit area is discretized into a grid with horizontal spacing of 1 meter and vertical spacing of 0.5 meters, with each grid node representing a target point. For each target point, the posterior probability density function is calculated using all hard and soft data according to the Bayesian conditional probability formula. The expected value of this posterior density function is the initial estimate of the parameter mean at that node. During the calculation, the uncertainty of the soft data is quantified by the variance of the regression residuals between it and the hard data, which is the average of the squared regression residuals at all borehole locations. This calculation is repeated for all grid nodes to obtain the initial distribution of the parameter mean across the entire pit area.

[0023] Then, the anisotropic variogram is calculated based on the initial distribution of parameter mean values, and the spatial autocorrelation distance is derived. The specific process is as follows: First, the values ​​of each grid node in the initial distribution of parameter mean values ​​are arranged according to spatial coordinates. For the horizontal direction, the experimental variogram is calculated along the length and width of the pit. Taking the length direction as an example, the node pairing distance is divided into several distance segments, each segment being 1 meter wide. The average of half the square of the parameter difference between all node pairs falling within the same distance segment is calculated to obtain the experimental variogram value corresponding to that distance segment. The experimental variograms for the width and vertical directions are calculated using the same method. Then, the experimental variogram is fitted using the spherical model theory. The expression of the spherical model is described as follows: when the distance is less than the range, the variogram equals the sill value multiplied by (1.5 times the distance divided by the range, then minus 0.5 times the cube of the distance divided by the cube of the range); when the distance is greater than or equal to the range, the variogram equals the sill value. The sill value is the maximum estimated value of the parameter variance, and the range is the spatial autocorrelation distance to be determined. By fitting the data using the least squares method, the sum of squared residuals between the theoretical and experimental variogram values ​​calculated by the spherical model is minimized, thereby determining the autocorrelation distance in the anisotropic direction. For the vertical direction, due to the limited number of data points, the maximum likelihood estimation method is used to directly invert the vertical autocorrelation distance.

[0024] Finally, the parameter variance is locally corrected using cross-validation residuals, and a prior random field is output. Specifically, one borehole location is reserved as a validation point. Using the hard data from the remaining boreholes and all soft data, the mean parameter value at this reserved location is re-estimated using the Bayesian maximum entropy method described above. The residual between the estimated and measured values ​​is calculated. This process is repeated for each borehole location to obtain a cross-validation residual sequence for all borehole locations. The variance of this residual sequence is calculated as the global residual variance. Then, for each grid node, a different correction coefficient is assigned based on its spatial distance from the nearest borehole: 0.8 for distances less than 5 meters, 1.0 for distances between 5 and 15 meters, and 1.2 for distances greater than 15 meters. The global residual variance is multiplied by the correction coefficient to obtain the local variance correction for that node. The initial parameter variance is then added to this local correction to obtain the corrected parameter variance. Thus, the mean of the parameters at each grid node (i.e. the mean in the initial distribution) and the corrected variance of the parameters, as well as the autocorrelation distances in the horizontal and vertical directions, are obtained. These three together constitute the prior random field.

[0025] In S2, Bayesian compressed sensing is used to perform sparse observation inversion on the prior random field, outputting the mean and covariance matrix of the parameter space under the posterior probability distribution, specifically including: First, the autocorrelation matrix of the prior random field is discretized into a sparse basis matrix to establish a linear mapping relationship between the observed data and the grid parameters. The mean parameter value at each grid node in the prior random field constitutes a high-dimensional vector, the dimension of which is equal to the total number of grid nodes (e.g., the number of nodes in the horizontal direction multiplied by the number of nodes in the vertical direction). Based on the autocorrelation distances and parameter variances obtained in the prior random field, an autocorrelation matrix is ​​constructed. Each element in the matrix corresponds to the autocorrelation value between any two grid nodes. This autocorrelation value is calculated using an exponential autocorrelation function: for any two nodes, their spatial distances in the horizontal and vertical directions are calculated, divided by the corresponding autocorrelation distances, and the sum of the squares is taken as the square root of the ratio. This ratio is used as the independent variable of the exponential autocorrelation function, and the function value is an exponent with the natural constant as the base and the independent variable being negative. The autocorrelation matrix is ​​decomposed using Cholleski to obtain a lower triangular matrix. Each column of this lower triangular matrix is ​​then normalized, and the eigenvectors corresponding to the K largest eigenvalues ​​are used to construct a sparse basis matrix. The value of K is determined based on the energy retention ratio, which is set to 95%, meaning the sum of the selected eigenvalues ​​accounts for more than 95% of the sum of all eigenvalues. The observation data includes estimated values ​​of sparse borehole geotechnical parameters and transformed continuous in-situ test data. These observation data are arranged into observation vectors. A linear relationship is established between the observation vectors and the grid parameter vectors through the sparse basis matrix: the grid parameter vector equals the sparse basis matrix multiplied by a sparse coefficient vector, and the observation vector equals the observation matrix multiplied by the grid parameter vector plus observation noise. Each row of the observation matrix corresponds to an observation point, and its elements are determined by the interpolation of the sparse basis matrix at the observation point location. The variance of the observation noise is determined based on the variance of the regression residuals between the continuous in-situ test data and the hard data.

[0026] Secondly, a Laplace prior constraint is introduced, and variational Bayesian inference is used to alternately update the hyperparameters and posterior covariance matrix. Specifically, the sparse coefficient vector is assigned a Laplace prior distribution, with its scale parameter initially set to one-tenth of the standard deviation of all observed data. The Laplace prior distribution can be decomposed into a Gaussian scale mixture form, where each sparse coefficient follows a normal distribution with a mean of zero and a variance equal to the square of an auxiliary scale parameter multiplied by a common variance parameter. This auxiliary scale parameter follows a gamma distribution with shape and scale parameters, where the shape parameter is set to 1 and the scale parameter is set to the square root of the Laplace prior distribution's scale parameter. Hyperparameters are defined as the common variance parameter and the observation noise variance, with their initial values ​​set to one-tenth and one-hundredth of the observed data variance, respectively. The variational Bayesian inference iterative update is employed, with the following steps: In the t-th iteration, the posterior distributions of the current hyperparameters and auxiliary scale parameters are fixed. The posterior covariance matrix is ​​calculated, which is equal to the transpose of the observation matrix multiplied by the observation matrix itself, multiplied by the reciprocal of the observation noise variance, plus a diagonal matrix. The diagonal elements of this diagonal matrix are equal to the reciprocals of the expected values ​​of each auxiliary scale parameter. The inverse of this sum matrix is ​​then used to obtain the posterior covariance matrix. Next, the posterior mean is calculated, which is equal to the posterior covariance matrix multiplied by the transpose of the observation matrix, multiplied by the observation vector, and divided by the observation noise variance. Then, the posterior expected values ​​of the auxiliary scale parameters are updated. The posterior expected value of the i-th auxiliary scale parameter is equal to the shape parameter plus one divided by the scale parameter plus the square of the i-th posterior mean divided by twice the common variance parameter, and then the reciprocal is taken. Finally, the hyperparameters are updated. The common variance parameter is updated to the arithmetic mean of the posterior expected values ​​of all auxiliary scale parameters, and the observation noise variance is updated to the observation vector minus the arithmetic mean of the observation matrix multiplied by the square of the posterior mean. The above iterative process is carried out sequentially. In each iteration, the posterior covariance matrix and posterior mean are updated first, then the posterior expectation of the auxiliary scale parameter is updated, and finally the hyperparameter is updated.

[0027] Finally, an iteration stopping condition is set: the relative change between the posterior means of two adjacent iterations is less than a set threshold. Specifically, the posterior mean vector obtained in the current iteration is subtracted from the posterior mean vector of the previous iteration, the L2 norm of the difference vector is calculated, and then divided by the L2 norm of the posterior mean vector of the previous iteration to obtain the relative change. The threshold is set to one-thousandth. Iteration stops when the relative change is less than one-thousandth, and the current posterior mean is output as the full grid mean of the parameter space, along with the posterior covariance matrix. The off-diagonal elements of this posterior covariance matrix reflect the correlation between parameters of different grid nodes, while the diagonal elements reflect the uncertainty variance of each node's parameters. This completes the Bayesian compressed sensing inversion, obtaining the full grid mean and covariance matrix of the parameter space under the posterior probability distribution.

[0028] In S3, multiple sets of random field implementations are generated based on the posterior probability distribution. Each set of random field implementations is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. Specifically, this includes: First, an incomplete Cholliski decomposition is performed on the covariance matrix under the posterior probability distribution to obtain a low-rank approximation lower triangular decomposition matrix. The posterior covariance matrix is ​​a square matrix whose dimension is equal to the total number of grid nodes after discretization of the pit area, denoted as N. The elements of this matrix have been obtained in step S2, and each element represents the posterior covariance of parameter values ​​at any two grid nodes. Since N is usually large (e.g., 50 nodes in the horizontal direction and 20 nodes in the vertical direction, then N=1000), the computational cost of directly performing Cholliski decomposition on the complete covariance matrix is ​​approximately N cubed, i.e., on the order of billions of times, which is computationally burdensome. Therefore, an incomplete Cholliski decomposition technique is adopted: a filling threshold is set, with a value of 0.01, and only lower triangular elements with absolute values ​​greater than this threshold are retained during the decomposition process, ignoring the remaining smaller elements. Specifically, the decomposition is performed row by row from row 1 to row N. For the current row, the square root of the diagonal elements of the row is calculated first, and then the decomposition value of each column (column number less than row number) in the row is calculated. The calculation formula is: ; in, The lower triangular decomposition matrix represents the first... Line number Column elements, Describe the posterior covariance matrix of the 1st generation. Line number Column elements, For the first The summation symbol represents the sum of elements from column 1 to column 2. The product of the elements decomposed from the column is accumulated. When the calculated... When the absolute value is less than the padding threshold of 0.01, it is directly set to zero and not stored. Through this process, a sparse lower triangular decomposition matrix is ​​obtained. Its number of non-zero elements is far less than the square of N, thus reducing the subsequent computation and storage requirements.

[0029] Secondly, conditional simulations are performed sequentially across all grid nodes to generate a set of random fields. This process utilizes the obtained full-grid mean vector of the parameter space (denoted as...). And the aforementioned lower triangular decomposition matrix L. First, generate an N-dimensional standard normal random vector, where each component is independent and follows a normal distribution with mean 0 and variance 1. Then, generate a temporary vector using the following formula: ; in, For temporary vectors, It is the mean vector of the entire grid (dimension N). The lower triangular matrix is ​​obtained from the incomplete Choreski decomposition. The above-generated standard normal random vector is denoted by " "" denotes matrix-vector multiplication. This formula transforms independent standard normal random vectors into multivariate normal random vectors with a specified posterior covariance structure. During the calculation, because... Since it is a sparse matrix, multiplication operations only involve non-zero elements, and the computational complexity is approximately the number of non-zero elements multiplied by N. The vector is essentially a random field implementation, where each component corresponds to a parameter value (e.g., elastic modulus) at a grid node. The process of generating a standard normal random vector and calculating is repeated. The process is repeated M times, with M being 2000 times, resulting in 2000 sets of random fields. No condition adjustment is needed during each generation process because the above formula implicitly includes the posterior covariance structure, thus achieving unconditional simulation. If conditional simulation is required (i.e., forcing the simulated value to equal the measured value at the borehole location), then after generating the standard normal random vector, the Kriging conditional simulation method is used for adjustment. Specifically, the difference between the simulated and measured values ​​at the borehole location is first calculated, and then this difference is propagated to all grid nodes through the covariance matrix for correction.

[0030] Next, for each set of random field implementations, corresponding excavation steps and support conditions were matched, and the finite element solver was called sequentially to calculate the spatial deformation field and time history curves of each construction stage. Specifically, the commercial finite element software PLAXIS 3D was used as the solver, automatically invoked through its Python script interface. First, an independent calculation directory was created for each set of random field implementations, named according to the implementation number from 1 to 2000. In each directory, based on the elastic modulus, cohesion, and other parameter values ​​provided by the random field implementation for each mesh node, linear interpolation was used to map them to each element integration point of the finite element model. The dimensions of the foundation pit geometry model were: width three times the foundation pit width, depth twice the foundation pit depth, and boundary conditions set as bottom fixed constraints and side normal constraints. The excavation steps were divided into six steps according to the actual construction plan, including: initial ground stress equilibrium, first support installation, first layer of earthwork excavation, second support installation, second layer of earthwork excavation, and pit bottom sealing. In this simulation, the support parameters (such as support stiffness and pile embedment depth) were initially set to a set of values ​​(e.g., pile diameter 1.0 meter, embedment depth 5 meters). The finite element solver was sequentially called to calculate the displacement field at the end of each step, recording the horizontal displacement and vertical settlement values ​​of each monitoring point (e.g., pit edge surface point, pile top, pile body), as well as the displacement field of all nodes as the spatial deformation field. Simultaneously, the construction time points corresponding to each step were recorded as time history curves. The numerical simulation of each step used a small-strain soil hardening constitutive model, and the convergence criterion was set to a displacement norm less than one-thousandth. After completing all step calculations for a set of random field implementations, the solver was closed, memory was cleared, and the calculation of the next implementation began. Due to the massive computational load of 2000 implementations, parallel computing was adopted, running simultaneously on 8 computing nodes, with each node processing one implementation at a time.

[0031] Finally, the grid node values ​​of the parameter space distribution are stored pairwise with the corresponding deformation field node values ​​to form a training sample set. For each random field implementation, its parameter space distribution is a vector of length N, recording the elastic modulus (or other geotechnical parameters) at each grid node. The corresponding deformation field output is the displacement value of all grid nodes after each monitoring step, which can also be organized into a vector of length N multiplied by the number of steps. The parameter vector is used as the input sample, and the deformation vector is used as the output label, stored in a file in binary format. The file header records the total number of samples (2000), the input dimension (N), and the output dimension (N multiplied by 6). In addition, the support parameters (initial values) and working condition identifier corresponding to each sample are also stored. The resulting training sample set contains 2000 input-output pairs, which are used in subsequent steps to construct the nonlinear mapping function.

[0032] In S4, a nonlinear mapping function from the parameter space distribution to the deformed field is constructed using the training sample set, and the prediction variance of the mapping output is quantified, specifically including: First, the parameter space distribution of the training sample set is used as input, and the deformation field is used as output to construct a hierarchical structure of a deep Gaussian process. Each layer employs an induced point sparsity approximation to reduce computational complexity. The deep Gaussian process consists of multiple stacked Gaussian process layers. The first layer uses the parameter space distribution vector as input, and its output serves as the input for the second layer, and so on, for a total of three layers. Each layer's Gaussian process uses a radial basis function kernel. The value of this kernel function is determined by exponentiating the negative square of the Euclidean distance between two input vectors divided by the length scale parameter, and then multiplying it by the signal variance parameter. To avoid directly calculating the kernel matrix between all training samples in each layer (with 2000 samples and a kernel matrix dimension of 2000 x 2000, the computational cost is on the order of millions), an induced point technique is introduced in each layer. Specifically, 200 points are uniformly selected from the input space of the training samples as induced points, and the initial values ​​of the induced point positions for each layer are obtained through random sampling. The kernel matrix of this layer is approximately represented as the product of the kernel matrix between induced points and the kernel matrix between induced points and training points, thus reducing the computational complexity from the cube of the number of samples to the square of the number of samples multiplied by the number of induced points.

[0033] Secondly, variational Bayesian inference is used to jointly learn the length scale, signal variance, and induced point positions of the kernel functions of each layer, while simultaneously optimizing the marginal likelihood variational lower bound. The training sample set is randomly divided into 10 subsets, each containing 200 training samples. The marginal likelihood variational lower bound is represented as the sum of the local lower bounds on each subset. The stochastic natural gradient descent algorithm is used to estimate the gradient of the lower bound with respect to the variational parameters. Specifically, in each iteration, a subset is randomly selected, the local lower bound corresponding to that subset is calculated, and then the natural gradient of the local lower bound with respect to the variational parameters (including the posterior mean and posterior covariance parameters of each Gaussian process) is calculated. The natural gradient is equal to the ordinary gradient multiplied by the inverse of the Fisher information matrix. The step size is set to 0.01, and the momentum coefficient is 0.9. After each update of the variational parameters, the length scale and signal variance of the current kernel function are fixed, and the induced point positions are updated using Riemann gradient ascent based on the current variational posterior covariance matrix. Riemann gradient ascent refers to calculating the gradient direction on the manifold space containing the induced point. This gradient is equal to the ordinary gradient of the objective function (i.e., the lower bound of the marginal likelihood variation) with respect to the coordinates of the induced point, minus the correction term along the tangent of the manifold. Specifically, the partial derivative of the lower bound with respect to the coordinates of each induced point is first calculated. This partial derivative is then projected onto the tangent space of the manifold to obtain the Riemann gradient. The coordinates of the induced point are then moved along the Riemann gradient direction with a step size of 0.005. The two updates are performed alternately: first, five consecutive iterations of variational parameter updates are performed, followed by one iteration of the induced point position update, and so on. The iteration stopping condition is set to the relative increment of the lower bound of the marginal likelihood variation between two adjacent iterations being less than a preset threshold, which is set to one-thousandth. The relative increment is calculated as: the lower bound value of the current iteration minus the lower bound value of the previous iteration, divided by the lower bound value of the previous iteration. The iteration stops when the relative increment is less than one-thousandth, and the kernel function parameters of each layer after joint learning (i.e., the length scale parameter and signal variance parameter of each layer) and the final spatial coordinates of 200 induced points in each layer are output.

[0034] Finally, based on the learned deep Gaussian process, the posterior prediction distribution is calculated for any parameter distribution in the input space, and the predicted mean and prediction covariance matrix of the output deformation field are obtained. The diagonal elements of the covariance matrix are extracted as the prediction variance. For a new input parameter distribution vector (denoted as the test input), it is first input into the first layer of the trained deep Gaussian process. This layer calculates the posterior output distribution corresponding to the test input based on the learned induced points, kernel function parameters, and variational posterior distribution. This distribution is a Gaussian distribution, and its mean and variance are given by the standard Gaussian process prediction formula. The specific calculation process is as follows: the kernel function values ​​between the test input and each induced point form a row vector, which is multiplied by the inverse of the kernel matrix between the induced points, and then multiplied by the posterior mean vector at the induced points to obtain the mean of the layer's output; the variance of the layer's output is equal to the self-kernel function value of the test input minus the above row vector multiplied by the inverse of the kernel matrix of the induced points and then multiplied by the transpose of the row vector, plus the contribution term of the variational posterior covariance. The output mean of this layer is used as the input to the next layer, and the above calculation is repeated until all three layers are passed. The final output mean of the third layer is the predicted mean vector of the deformed field (dimension equal to the number of output nodes), and the output variance vector of the third layer is the predicted variance at each node of the deformed field. The predicted variance is used as the diagonal elements, while the off-diagonal elements are not output for the time being, forming the predicted variance vector. This completes the construction of the nonlinear mapping from the parameter space distribution to the deformed field and the quantization of the predicted variance.

[0035] In S5, using deformation control threshold and allowable exceedance probability as constraints, a reverse search is performed on the support parameter combination to ensure that the probability of satisfying the constraints under the deformation probability distribution predicted by the nonlinear mapping function is higher than a preset level. Specifically, this includes: First, based on the deformation mean and prediction variance output by the nonlinear mapping function, an expected improvement acquisition function is constructed, which includes deformation control threshold constraints and over-limit penalties. The deformation control threshold is determined according to the surrounding environment of the foundation pit; for example, when near a subway tunnel, the horizontal displacement threshold is set to 5 mm. For each set of candidate support parameters (including pile diameter, embedment depth, anchor cable length, and support stiffness), the predicted mean and prediction variance of each control point in the deformation field are calculated using the depth Gaussian process in step S4. The over-limit probability is defined as the probability that the deformation is greater than the control threshold. This probability is calculated using the standard normal cumulative distribution function: subtract the predicted mean from the control threshold, then divide by the predicted standard deviation (i.e., the square root of the predicted variance) to obtain the standardized deviation. Then, the over-limit probability value is obtained by looking up a table using the standard normal cumulative distribution function. The expected improvement acquisition function consists of two parts: the first part is the expected improvement of the deformation margin under the current optimal support cost, and the second part is the over-limit penalty term. The over-limit penalty term is equal to the difference between the over-limit probability and the preset level (e.g., 5%) multiplied by the penalty coefficient. When the over-limit probability is less than the preset level, the penalty term is zero. A larger acquisition function value indicates a better combination of candidate support parameters.

[0036] Secondly, an adaptive evolutionary strategy based on the covariance matrix is ​​employed for iterative searching within the support parameter space. The ranges of each dimension of the support parameter space are set based on engineering experience: pile diameter ranges from 0.6 meters to 1.5 meters, embedment depth ranges from 3 meters to 10 meters, anchor cable length ranges from 10 meters to 25 meters, and support stiffness ranges from 10 meganewtons per meter to 100 meganewtons per meter. The initial population size of the evolutionary strategy is 16 candidate solutions, each candidate solution being a four-dimensional vector (corresponding to the four support parameters mentioned above). The initial covariance matrix is ​​set as an identity matrix, and the initial step size is one-tenth of the width of each dimension. In each iteration, a new population of candidate solutions is generated based on the current mean and covariance matrix. For each candidate solution, the depth Gaussian process in step S4 is called to calculate the deformation mean and prediction variance, and then the out-of-limit probability and the expected improved acquisition function value are calculated. Based on the acquisition function values, the top 8 optimal candidate solutions are selected to update the mean vector, covariance matrix, and step size. During the update process, the penalty weight is dynamically adjusted based on the difference between the out-of-limit probability of each candidate solution and a preset level (5%): when the out-of-limit probability is greater than the preset level, the acquisition function value of the candidate solution is multiplied by a decay factor, with an initial value of 0.8, which is linearly reduced to 0.2 according to the amount of out-of-limit probability; when the out-of-limit probability is less than or equal to the preset level, the penalty weight is set to 1 and does not decay. The maximum value of the acquisition function in each iteration is also recorded.

[0037] Finally, the iteration termination condition is set as follows: the relative change of the optimal value of the acquisition function in 20 consecutive iterations is less than a set tolerance. The relative change is calculated as follows: the optimal value of the current iteration minus the optimal value 20 iterations ago, then divided by the optimal value 20 iterations ago. The tolerance is set to 0.1%. The search terminates when the absolute value of the relative change in 20 consecutive iterations is less than 0.1%. From the candidate solutions of the last iteration, candidate solutions with an out-of-limit probability of less than or equal to 5% are selected, and the one with the lowest support cost is chosen (the support cost is estimated by linear weighting of the product of pile length and pile diameter, anchor cable length, and support stiffness, with each weight coefficient determined by engineering economic indicators). The support parameter combination corresponding to this candidate solution is output as the final optimization result.

[0038] Please see Figure 2 As shown, a system for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters includes: The data acquisition and prior modeling module collects sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area, and constructs a prior random field including parameter mean, variance and autocorrelation distance. The Bayesian compressed sensing inversion module uses Bayesian compressed sensing to perform sparse observation inversion on a prior random field and outputs the mean and covariance matrix of the parameter space under the posterior probability distribution. The random field generation and numerical simulation module generates multiple sets of random field implementations based on the posterior probability distribution. Each set of random field implementations is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. The nonlinear mapping learning module uses the training sample set to construct a nonlinear mapping function from the parameter space distribution to the deformation field, and quantifies the prediction variance of the mapping output. The robust inverse optimization module uses deformation control threshold and allowable over-limit probability as constraints to inversely search for support parameter combinations, so that the probability of satisfying the constraints under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

[0039] The working principle of this invention is as follows: First, sparse borehole geotechnical parameters and continuous in-situ test data of the foundation pit area are collected to construct a prior random field containing parameter mean, variance, and autocorrelation distance. Then, Bayesian compressed sensing is used to perform sparse observation inversion on the prior random field, outputting the full-grid mean and covariance matrix of the parameter space under the posterior probability distribution. Next, multiple sets of random fields are generated based on the posterior probability distribution, and they are successively substituted into finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter space distribution and the deformation field. Then, the training sample set is used to construct a nonlinear mapping function from the parameter space distribution to the deformation field, and the predicted variance of the mapping output is quantified. Finally, with deformation control threshold and allowable exceedance probability as constraints, the combination of support parameters is searched in reverse so that the probability of satisfying the constraint conditions under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

[0040] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters, characterized in that, Includes the following steps: S1. Collect sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area to construct a prior random field containing parameter mean, variance and autocorrelation distance; S2 uses Bayesian compressed sensing to perform sparse observation inversion on the prior random field and outputs the mean and covariance matrix of the parameter space full grid under the posterior probability distribution. S3, based on the posterior probability distribution, generates multiple sets of random fields. Each set of random fields is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. S4. Construct a nonlinear mapping function from the parameter space distribution to the deformable field using the training sample set, and quantify the prediction variance of the mapping output. S5 uses the deformation control threshold and the allowable over-limit probability as constraints to inversely search for the support parameter combination, so that the probability of satisfying the constraint conditions under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

2. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 1, characterized in that, S1 specifically includes: The sparse borehole geotechnical parameters are used as hard data, and the continuous in-situ test data are used as soft data. The initial distribution of the parameter mean is generated by fusion. The anisotropic variability function is calculated based on the parameter mean distribution, and the spatial autocorrelation distance is then derived. The parameter variance is locally corrected using cross-validation residuals, and a prior random field is output.

3. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 1, characterized in that, S2 specifically includes: Discretize the autocorrelation matrix of the prior random field into a sparse basis matrix and establish a linear mapping relationship between the observed data and the grid parameters; We introduce Laplace prior constraints and use variational Bayesian inference to alternately update hyperparameters and posterior covariance matrices. Stop when the change in the posterior mean between two consecutive iterations is less than a set threshold, and output the full grid mean and covariance matrix of the parameter space under the posterior probability distribution.

4. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 3, characterized in that, The introduction of Laplace prior constraints and the use of variational Bayesian inference to alternately update hyperparameters and posterior covariance matrices specifically include: A sparse linear observation equation is constructed based on the autocorrelation matrix of the prior random field and the observation data, and a Laplace prior distribution is assigned to the full grid mean of the parameter space. The Laplace prior distribution is decomposed into a Gaussian scale mixture, and an auxiliary scale parameter is introduced and its gamma prior is set. The variational Bayesian expectation-maximization step is adopted to sequentially update the posterior expectation, posterior covariance matrix and hyperparameter of the auxiliary scaling parameter until the difference between the posterior means of two adjacent iterations is less than the preset tolerance, and then the updated posterior mean and covariance matrix are output.

5. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 1, characterized in that, S3 specifically includes: The covariance matrix under the posterior probability distribution is decomposed by Cholliski decomposition, and combined with the full grid mean, a predetermined number of random fields are generated by sequential Gaussian simulation. For each set of random fields, the corresponding excavation sequence and support conditions are matched, and the finite element solver is called sequentially to calculate the spatial deformation field and time history curve of each construction stage. The grid node values ​​of the parameter space distribution are stored one pair at a time with the corresponding deformation field node values ​​to form a training sample set.

6. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 5, characterized in that, The generation of a predetermined number of random fields is specifically implemented as follows: The posterior covariance matrix is ​​subjected to incomplete Choleski decomposition to obtain a low-rank approximation lower triangular decomposition matrix. Conditional simulations are performed sequentially across all grid nodes. After each node's parameter value is generated, the covariance conditions of adjacent unsampled nodes are immediately adjusted based on the updated values ​​of the corresponding nodes. Repeat the sequential process until all grid nodes are assigned values, and output a set of random field implementations that follow a posterior probability distribution.

7. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 1, characterized in that, S4 specifically includes: The parameter spatial distribution in the training sample set is used as input and the deformation field is used as output to construct a hierarchical structure of the deep Gaussian process. Each layer adopts the induced point sparsity approximation to reduce the computational complexity. The length scale, signal variance, and induced point location of the kernel function of each layer are jointly learned by variational Bayesian inference, while the variational lower bound of the marginal likelihood is optimized. Based on the learned deep Gaussian process, the posterior prediction distribution is calculated for any parameter distribution in the input space, and the predicted mean and prediction covariance matrix of the output deformation field are obtained. The diagonal elements of the covariance matrix are extracted as the prediction variance.

8. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 7, characterized in that, The method of inferring the length scale, signal variance, and induced point location of the kernel function of each layer through variational Bayesian inference, while simultaneously optimizing the variational lower bound of the marginal likelihood, specifically includes: The marginal likelihood variational lower bound is decomposed into the sum of local lower bounds on each training subset, and the gradient of the lower bound with respect to the variational parameters is estimated by stochastic natural gradient descent. With the kernel function's length scale and signal variance fixed, the Riemann gradient ascent update is performed on the induced point position based on the current variational posterior covariance matrix. Alternately update the kernel function hyperparameters and induced point positions until the relative increment of the variational lower bound between two adjacent iterations is lower than a preset threshold, and output the kernel function parameters and induced point positions of each layer after joint learning.

9. The method for predicting the spatiotemporal evolution of foundation pit deformation and for inverse optimization of support parameters according to claim 1, characterized in that, S5 specifically includes: Based on the deformation mean and prediction variance output by the nonlinear mapping function, an expected improved acquisition function is constructed, which includes deformation control threshold constraints and over-limit penalties. An adaptive evolution strategy based on the covariance matrix is ​​used to iteratively search within the support parameter space. In each iteration, the acquisition function is called to evaluate the deformation over-limit probability corresponding to the candidate parameter combination, and the penalty weight is dynamically adjusted based on the difference between the over-limit probability and the preset level. The search terminates when the relative change of the optimal value of the acquisition function is less than the set tolerance in multiple consecutive iterations, and outputs the combination of support parameters that meets the over-limit probability requirement and minimizes the support cost.

10. A system for predicting the spatiotemporal evolution of foundation pit deformation and for inversely optimizing support parameters, characterized in that, A method for performing the spatiotemporal evolution prediction of foundation pit deformation and inverse optimization of support parameters as described in any one of claims 1-9 includes: The data acquisition and prior modeling module collects sparse borehole geotechnical parameters and continuous in-situ test data in the foundation pit area, and constructs a prior random field including parameter mean, variance and autocorrelation distance. The Bayesian compressed sensing inversion module uses Bayesian compressed sensing to perform sparse observation inversion on a prior random field and outputs the mean and covariance matrix of the parameter space under the posterior probability distribution. The random field generation and numerical simulation module generates multiple sets of random field implementations based on the posterior probability distribution. Each set of random field implementations is then substituted into the finite element numerical simulation to calculate the corresponding spatiotemporal deformation response of the foundation pit, forming a training sample set between the parameter spatial distribution and the deformation field. The nonlinear mapping learning module uses the training sample set to construct a nonlinear mapping function from the parameter space distribution to the deformation field, and quantifies the prediction variance of the mapping output. The robust inverse optimization module uses deformation control threshold and allowable over-limit probability as constraints to inversely search for support parameter combinations, so that the probability of satisfying the constraints under the deformation probability distribution predicted by the nonlinear mapping function is higher than the preset level.

Images