A random field construction and engineering analysis method considering uncertainty of geotechnical parameters

Abstract

This invention discloses a method for constructing and analyzing random fields considering the uncertainties of soil and rock parameters, belonging to the field of numerical analysis in geotechnical engineering. This invention hierarchically organizes the three-dimensional coordinates and physical and mechanical parameters of the site, determines the optimal marginal distribution of key soil parameters through Bayesian inference, constructs a deterministic finite element numerical model, extracts the three-dimensional center point coordinates of the grid cells, fits the autocorrelation function and Copula model to obtain spatial correlation distance and parameter correlation coefficients, constructs a joint covariance matrix, and outputs a three-dimensional conditional random field through Karhunen-Loève expansion and joint ordinary kriging condition correction. This invention can effectively quantify the statistical uncertainty of soil parameters, overcome the problem of insufficient accuracy in parameter distribution identification under small sample exploration conditions, characterize the spatial anisotropy of soil and rock parameters and the nonlinear correlation of multiple parameters, and significantly improve the realism of three-dimensional conditional random field modeling.

Description

Technical Field

[0001] This invention relates to the field of numerical analysis technology in geotechnical engineering, and specifically to a method for constructing and analyzing random fields that considers the uncertainty of geotechnical parameters. Background Technology

[0002] Soil and rock masses are naturally formed heterogeneous geological bodies, and their physical and mechanical parameters exhibit significant spatial variability and uncertainty, which are core factors affecting the reliability of numerical analysis results in geotechnical engineering. Random field theory is currently the mainstream technique for characterizing the spatial variability of soil and rock parameters. By modeling soil and rock parameters as spatial random functions, the spatial distribution law of the parameters can be quantitatively described, providing data support for reliability analysis and risk assessment in geotechnical engineering.

[0003] Publication No. CN115357994A discloses a spatial random field modeling method for surrounding rock parameters. This method selects only two types of deformation parameters, elastic modulus and Poisson, to construct a two-dimensional joint distribution. It also relies solely on the KS test criterion to identify the edge distribution, which cannot quantify the statistical uncertainty of the parameters. This makes it difficult to adapt to complex geotechnical engineering multi-parameter coupled analysis scenarios. Under small sample exploration conditions, the accuracy of parameter distribution identification is low, which greatly reduces the accuracy and reliability of random field modeling. Summary of the Invention

[0004] In view of this, the present invention provides a method for constructing and analyzing random fields considering the uncertainty of geotechnical parameters, integrating Bayesian MCMC parameter inference, multi-dimensional correlation structure fitting, and measured constraint 3D conditional random field construction into a unified process. This invention leverages the small-sample identification advantage of Bayesian methods, combines on-site measured data to complete parameter correction, quantifies the uncertainty of geotechnical parameters, and effectively improves the accuracy of 3D conditional random field modeling.

[0005] The technical solution of this invention is implemented as follows: This invention provides a method for constructing and analyzing random fields considering the uncertainty of geotechnical parameters, including the following steps: S1: Collect three-dimensional coordinate data and physical and mechanical parameter data of each soil layer in the site, divide the soil layers according to depth and collect them in layers, calculate the statistical characteristics of the parameters of each layer, determine the key soil parameters, and determine the optimal marginal distribution of each key soil parameter and the posterior estimate of its distribution parameters based on Bayesian inference. S2: Construct a deterministic finite element numerical model using the average values ​​of physical and mechanical parameters, extract the element topology information and three-dimensional coordinate information of the mesh nodes of the model, and calculate the three-dimensional center point coordinates of each mesh element; S3: Construct a depth data sequence using the three-dimensional coordinate data from step S1, calculate the empirical autocorrelation value, fit the autocorrelation function to determine the spatial correlation distance in each direction, and combine the optimal edge distribution determined in step S1 to fit the Copula model to obtain the correlation coefficient between parameters. S4: Based on the coordinates of the three-dimensional center point of the grid cell, the spatial correlation distance, the correlation coefficient of the optimal Copula model, and the distribution parameters estimated from the posterior, construct the joint covariance matrix of the key soil parameters, use Karhunen-Loève expansion to generate unconditional random samples, combine the physical and mechanical parameter data from step S1 to perform condition correction, and output the three-dimensional conditional random field of the key soil parameters.

[0006] Based on the above technical solutions, preferably, the physical and mechanical parameter data in step S1 are obtained through drilling sampling, indoor geotechnical tests, in-situ detection and wave velocity testing. The in-situ detection and wave velocity test data both have continuous depth coordinates. The statistical characteristics include maximum value, minimum value, average value, standard deviation and coefficient of variation. The key soil parameters are the parameters with larger coefficients of variation among the physical and mechanical parameters of each soil layer.

[0007] Based on the above technical solutions, the preferred procedure for Bayesian inference and posterior estimation in step S1 is as follows: four candidate probability distribution models are selected, and a uniform prior distribution is set for the distribution parameters of each model; the likelihood function of each model is constructed based on the data of key soil parameters; according to Bayes' theorem, the posterior distribution of the distribution parameters is calculated from the prior distribution and the likelihood function; the model evidence of each model is calculated by Monte Carlo integration, and the posterior probability of each model is calculated based on the model evidence, and the model with the largest posterior probability is selected as the optimal marginal distribution; the distribution parameters of the optimal marginal distribution are posteriorly sampled using the Metropolis-Hastings Markov chain Monte Carlo method to obtain the posterior estimate of the distribution parameters, and a prediction sample is generated through the posterior estimate, and the fit of the optimal marginal distribution is verified by comparing it with the measured data.

[0008] More preferably, the alternative probability distribution models are normal distribution, log-normal distribution, gamma distribution, and Weibull distribution. All four types of distributions include location parameters and scale parameters, with gamma and Weibull distributions additionally including shape parameters.

[0009] More preferably, the uniform prior interval of the position parameter is [0, 2]. m The uniform prior interval for scale and shape parameters is [0.1]. s ,2 s ], m This is the average value. s To eliminate the standard deviation, when calculating model evidence, a logarithmic transformation is used to subtract the maximum value to eliminate the underflow of Monte Carlo integral values.

[0010] Based on the above technical solutions, preferably, in step S2, when constructing the deterministic finite element numerical model, the average value of the physical and mechanical parameters is used as the fixed parameter, and boundary conditions and loads are set. The element topology information and the three-dimensional coordinate information of the grid nodes of the model are extracted. The element topology information is the grid mapping data of the connection relationship between elements and nodes. The associated nodes corresponding to each grid element are retrieved using the element topology information, and the three-dimensional coordinate information of the grid nodes corresponding to the associated nodes is read and the average value is calculated. The three-dimensional center point coordinates of the grid element are thus calculated.

[0011] Based on the above technical solutions, preferably, the specific process of autocorrelation function fitting and Copula model fitting in step S3 is as follows: determine the equal-interval sampling step size, calculate the empirical autocorrelation function values ​​of each physical and mechanical parameter under different spatial lag distances; use multiple autocorrelation function models for fitting, select the optimal autocorrelation function model based on the goodness of fit, calculate the spatial correlation distance of each soil layer's physical and mechanical parameters in the vertical and horizontal directions; use multiple Copula models to fit the correlation structure between each key soil parameter, select the optimal Copula model based on the fitting index, and output its correlation coefficient.

[0012] More preferably, the autocorrelation function model uses three types: squared exponential, single exponential, and cosine exponential. The positive values ​​of the empirical autocorrelation function are fitted using the Levenberg-Marquardt nonlinear least squares algorithm. When the goodness of fit of the three models is... R 2 When the value is less than 0.1, the relevant distances in the vertical and horizontal directions are determined by engineering experience.

[0013] More preferably, the process of fitting the optimal Copula model involves: obtaining the cumulative distribution function corresponding to the optimal edge distribution based on the optimal edge distribution and its distribution parameters determined in step S1, and mapping the data of each key soil parameter to a uniform space through the cumulative distribution function; adjusting the boundaries of the mapped data to ensure it is located within […]. e ,1- e Within the interval, e =1×10 -8 The correlation structure between key soil parameters was fitted using four types of Copula models: Gaussian, Frank, Gumbel, and t-Copula. The AIC (Akaike Information Criterion) was used to validate the models, and the Copula model with the smallest AIC value was selected as the optimal Copula model. Its correlation coefficient was then output.

[0014] Based on the above technical solutions, the preferred step S4 is as follows: An autocovariance matrix is ​​constructed based on the coordinates of the center point of the grid cell and the spatial correlation distances in each direction. The amplitude is determined by the posterior-estimated distribution parameters. A cross-covariance matrix between parameters is constructed using the optimal Copula correlation coefficient, forming a joint covariance matrix of the key soil parameters. Its eigenvalues ​​are decomposed, truncated to a cumulative characteristic energy ratio of not less than 0.95 while retaining the principal characteristic modes. An unconditional random sample is generated using Karhunen-Loève expansion. Conditional corrections are performed using the joint ordinary kriging method and physical and mechanical parameter data, outputting a three-dimensional conditional random field of the key soil parameters.

[0015] The beneficial effects of this invention are as follows: (1) This invention establishes a complete Bayesian posterior distribution inference system. It sets a layered uniform prior for small sample soil and rock data and uses Monte Carlo integral to calculate model evidence to screen the optimal marginal distribution. It can fully explore the statistical regularity in the limited exploration data, quantify the fluctuation of the parameters themselves, make the marginal distribution judgment more accurate, and provide reliable parameter model support for the construction of three-dimensional conditional random fields.

[0016] (2) This invention combines multiple autocorrelation functions with multiple types of Copula models for fitting and screening, which can comprehensively reflect the anisotropic characteristics of soil parameters such as cohesion, internal friction angle, and compression modulus in three-dimensional space, as well as the nonlinear correlation between parameters, and accurately calculate the spatial correlation distance in each direction and the correlation coefficient matrix between parameters, providing a dual basis for the construction of three-dimensional conditional random fields in terms of spatial variability and parameter synergy.

[0017] (3) The present invention first generates a random field through Karhunen-Loève expansion, and then uses borehole measured data to perform joint ordinary kriging condition correction to establish a three-dimensional conditional random field with measured data constraints, so that the parameter assignment of each grid unit matches the actual exploration information, better restores the natural spatial variation law of the rock and soil, and makes the three-dimensional conditional random field model closer to the actual situation on site. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a complete flowchart of the overall implementation of the method of the present invention; Figure 2This is a comparison chart of the posterior probabilities of the cohesion, internal friction angle, and compression modulus of the 3-2 layer silty clay of the present invention under four alternative probability distribution models. Figure 3 This is a schematic diagram of the deterministic finite element foundation model of the Wuhan rail transit tunnel section according to the present invention; Figure 4 This is a schematic diagram showing the variation of the horizontal stress index with depth obtained from the site flat shovel lateral expansion test of the present invention; Figure 5 The figure shows a comparison of the correlation fitting effects between cohesion and internal friction angle based on Gaussian, Frank, Gumbel, and t-Copula functions in this invention. (a) to (d) are the fitting results of the four types of Copula models, respectively. Figure 6 This is a spatial distribution cloud map of the compressibility modulus of a three-dimensional conditional random field generated by the Karhunen-Loève expansion and joint ordinary kriging correction according to the present invention. Detailed Implementation

[0020] To better understand the purpose, technical solutions, and advantages of this application, the application has been described and illustrated below with reference to the accompanying drawings and embodiments. However, those skilled in the art should understand that this application can be implemented without these details. It will be apparent to those skilled in the art that various modifications can be made to the embodiments disclosed in this application, and the general principles defined in this application can be applied to other embodiments and application scenarios without departing from the principles and scope of this application. Therefore, this application is not limited to the illustrated embodiments, but conforms to the broadest scope consistent with the scope of protection claimed in this application.

[0021] This embodiment uses the shield tunnel project between Science Park Station and Gangdu Garden Station on Wuhan Metro Line 12 as an application scenario to fully illustrate the method of the present invention. The overall process of the method is as follows: Figure 1 As shown.

[0022] S1: Parametric Statistics and Bayesian Distribution Identification 1.1 Parameter Collection and Screening of Key Soil Parameters The physical and mechanical parameters of the site soil layers were obtained through drilling sampling, laboratory geotechnical analysis, in-situ exploration, and wave velocity testing. The in-situ exploration and wave velocity test data included continuous depth coordinates. After dividing the soil into layers according to depth, the measured data of the two types of parameters were collected layer by layer, and the maximum, minimum, and average values ​​were calculated for each layer. m Standard deviation s Coefficient of variation (COV); the coefficient of variation is the dimensionless ratio of the standard deviation to the mean (COV = s / mThe main soil layers in this embodiment include miscellaneous fill, silty clay, and silty clay. The statistical results of the physical parameters of each soil layer are shown in Table 1, and the statistical results of the mechanical parameters of each soil layer are shown in Table 2.

[0023] Table 1. Statistics of physical parameters of soil layers in the site.

[0024] Table 2. Statistics of mechanical parameters of soil layers at the site

[0025] Comparing the coefficients of variation values ​​of all soil layers in Tables 1 and 2, the overall dispersion of mechanical parameters is higher than that of physical parameters. Among them, the coefficients of variation of three mechanical parameters—cohesion, internal friction angle, and compression modulus—of the silty clay in layer 3-2 are the largest. These three mechanical parameters are identified as the key soil parameters in this embodiment. Due to the limited disturbance range of the shield tunnel, only four representative soil layers within the project area were selected, and the coefficients of variation of the key soil parameters of each layer were statistically analyzed separately. The statistical results are summarized in Table 3.

[0026] Table 3. Statistics of coefficients of variation of key soil parameters for each soil layer

[0027] Table 3 shows that the coefficients of variation for all key soil parameters in layer 3-2 silty clay are generally high. Combined with Table 2, which indicates a sufficient number of experimental samples and stable statistical patterns, layer 3-2 silty clay was selected as the core stratum for this random field modeling. The measured key soil parameters of this stratum were collected to form a sample set Data, with each sample being Data. i ( i =1, 2, ..., N sam ), N sam This represents the total sample size.

[0028] 1.2 Bayesian Optimal Marginal Distribution Identification and Posterior Sampling 1.2.1 Selection of Alternative Distributions and Uniform Prior Setting Four candidate probability distribution models were selected, and a uniform prior distribution was set for the distribution parameters of each candidate probability distribution model. The reasonable range of values ​​for the distribution parameters was defined, providing the parameter integral domain for Bayesian inference of key soil parameters.

[0029] Four commonly used geotechnical engineering models—normal distribution, log-normal distribution, gamma distribution, and Weibull distribution—were selected as candidate probability distribution models, and a uniform prior distribution was assigned to the distribution parameters of each model. The distribution parameter vectors corresponding to each candidate probability distribution model are shown below. i The definition is as follows: normal distribution: i =[ m , s ]; Log-normal distribution: i =[ m ln , s ln ],in m ln The average of the parameters after taking the natural logarithm. s ln The standard deviation of the parameter after taking the natural logarithm; Gamma distribution: i =[ α , β ],in α For shape parameters, β For scale parameters; Weibull distribution: i =[ k , l ],in k For shape parameters, l This is the scale parameter.

[0030] Due to the limited sample size of the site soil mechanics test and the lack of sufficient prior experience to constrain the distribution parameters, this embodiment uniformly adopts a uniform uninformed prior for all distribution parameters and formulates parameter value boundary rules: the value range of the average location parameter is set to [0, 2]. m The value range for scale and shape parameters is set to [0.1]. s ,2 s ].average value m and standard deviation s The key soil parameters were uniformly taken from the sample statistics of the silty clay layer 3-2 in Table 2. Specific prior intervals for each parameter: cohesion. m ∈[0,34], s ∈[0.357,7.14]; Angle of internal friction m ∈[0,18], s ∈[0.162,3.24]; compressibility modulus m ∈[0,8.2], s ∈[0.043,0.86].

[0031] 1.2.2 Constructing the likelihood function of alternative models Based on the sample data of key soil parameters, the likelihood function corresponding to each alternative probability distribution model is constructed. This function is used to characterize the probability of the measured sample appearing under the given alternative probability distribution model and distribution parameters. It is the basic function for carrying out Bayesian inference.

[0032] This embodiment relies on the measured sample set Data of key soil parameters to construct the sample joint likelihood function corresponding to each alternative probability distribution model. For the first... i Single sample data i In the given first n Alternative probability distribution model M n With distribution parameter vector i When using the alternative probability distribution model, the probability density corresponding to that model is adopted. Characterizing the single-sample likelihood. Assuming all single samples are independent, the joint likelihood function is constructed according to the joint probability product rule of independent random variables. Its expression is the product of the single-sample probability densities:

[0033] In the formula: P (Data| i , M n ) represents the joint likelihood value.

[0034] 1.2.3 Calculating the posterior distribution of parameters based on Bayes' theorem When the number of physical and mechanical parameters is small, traditional fitting tests show poor stability in selecting distribution models. This invention, based on Bayes' theorem and combining a uniform prior distribution with the likelihood function, calculates the posterior distribution of the distribution parameters. Since the distribution parameters are continuous variables, discrete summation cannot be used to calculate the total probability. Instead, global integration eliminates the dependence on the distribution parameters, calculates model evidence, and uses this model evidence to quantitatively select the best among multiple alternative probability distribution models.

[0035] In this embodiment, the total number of probability distribution models available for screening and fitting key soil parameters is... N m =4, and any alternative probability distribution model is denoted as . M n ( n =1, 2... N m By Bayes' theorem, given a sample set of measured key soil parameters (Data), the posterior probability expression for the validity of this alternative probability distribution model is:

[0036] In the formula: P (Data |M n () serves as model evidence, used to measure the degree of fit between the model and the measured sample data; P ( M n) represents the prior probability of the model. When there is no prior experience, each model has equal probability. The denominator is a normalization constant, which will not change the relative ranking of the different models.

[0037] Total number of candidate models in this embodiment N m =4, each model is assigned an equal prior probability, that is P ( M 1)= P ( M 2)= P ( M 3)= P ( M 4) = 1 / 4, substituting the prior conditions of the equal model into the posterior probability expression, the constant terms in the numerator and denominator... They can cancel each other out, resulting in a simplified expression under the prior assumptions of the equal model:

[0038] It can be deduced that, under the premise that the candidate probability distribution models are a priori equally probable, P ( M n | Data) and P (Data |M n The optimal distribution model is the alternative probability distribution model that is directly proportional to the model evidence and achieves the maximum value.

[0039] The essence of model evidence is the marginal probability of a sample obtained after marginalizing the distribution parameters in the parameter space. If the distribution parameters are discrete variables, the total probability formula in a summation form can be used; however, in this example, the distribution parameters such as the mean and standard deviation are continuous variables, so discrete summation is no longer applicable. The total probability formula needs to be transformed into a global integral form within the parameter domain, i.e., the calculation formula for model evidence:

[0040] In the formula: Ω is the domain of the integral of the effective values ​​of the distributed parameter; P ( θ|M n ) represents the prior density of the distribution parameter.

[0041] When the sample data of physical and mechanical parameters is finite, a uniform, uninformative prior is used. Taking the two-parameter normal distribution and the log-normal distribution as examples, the prior density of the distribution parameters of the two is uniformly expressed as:

[0042] In the formula: i =[ m, s ], L b ,U b Let be the lower and upper bounds of the unified integral of the distributed parameter vector, respectively. L b =[ m min , s min ], U b =[ m max , s max ], m min , s min , m max , s max Values ​​are taken according to the prior interval rules in step 1.2.1. When m ∈[ m min , m max ]and s ∈[ s min , s max When ], the distribution of the parameter vector to be estimated i fall into[ L b , U b ], the prior density is 1 / ( U b - L b ); otherwise, take 0.

[0043] Based on the measured statistical values ​​of silty clay in the 3-2 layers obtained in step S1, the specific range of distribution parameters for each alternative probability distribution model is determined, thereby constructing a uniform prior density and conducting global integral calculation model evidence.

[0044] 1.2.4 Determining the Optimal Marginal Distribution through Monte Carlo Integration Since no analytical solution exists for the model evidence, the Monte Carlo numerical integration method is used to calculate the model evidence for each candidate probability distribution model. Then, based on the model evidence, the posterior probability of each candidate probability distribution model is calculated, and the candidate probability distribution model with the highest posterior probability is selected as the optimal marginal distribution. This embodiment extracts values ​​based on the distribution parameter ranges of each candidate probability distribution model set in step 1.2.1. N mc =100,000 sets of distributed parameter vectors i k ,in N mcThe total number of Monte Carlo integral samples; the distribution parameter vector for each group. i k Substituting the sample joint likelihood function expression constructed in step 1.2.2, we obtain the joint likelihood value corresponding to this set of parameters, and then take its natural logarithm to obtain the log-likelihood value ln. L k .

[0045] Directly using the summation of joint likelihood values ​​can easily produce extremely small floating-point numbers, leading to numerical underflow. To eliminate Monte Carlo integral numerical underflow, a method of subtracting the maximum value using a logarithmic transformation is employed. The calculation formula is as follows:

[0046] Where: max k (ln L k ) represents the maximum log-likelihood value corresponding to the distribution parameter vector obtained by sampling from the range of values ​​of the self-distributed parameters in each group; exp(·) is the natural exponential function, and exp(·) is uniformly defined throughout the text. x )=e x .

[0047] This method eliminates numerical underflow and improves the stability of model evidence calculation under small sample conditions. Using the calculated model evidence corresponding to each alternative distribution model, and combining it with Bayesian theory, the normalized posterior probability transformation is completed.

[0048] The candidate probability distribution model with the highest posterior probability was selected as the optimal marginal distribution for key soil parameters. In this embodiment, measured statistical data from 3-2 layers of silty clay were used as input. The optimal marginal distribution for cohesion (mean 17 kPa, standard deviation 3.57 kPa) was a log-normal distribution; the optimal marginal distribution for internal friction angle (mean 9°, standard deviation 1.62°) was a normal distribution; and the optimal marginal distribution for compression modulus (mean 4.1 MPa, standard deviation 0.43 MPa) was a Weibull distribution. The fitting results for each parameter are shown in […]. Figure 2 The optimal distribution type and distribution parameters are output and transmitted to steps S3 and S4 as calculation inputs.

[0049] 1.2.5 Posterior Sampling and Estimation of Distribution Parameters for the Optimal Marginal Distribution After determining the optimal marginal distribution type, posterior sampling of the distribution parameters of the optimal marginal distribution is required to obtain the posterior estimation results of the distribution parameters. This embodiment uses the Metropolis-Hastings Markov Chain Monte Carlo (MH-MCMC) algorithm to perform posterior sampling of the distribution parameter vector of the optimal marginal distribution. The candidate parameter generation formula is as follows:

[0050] In the formula: i prop This is a candidate parameter vector; i curr This is the current iteration parameter vector; s prop The suggested standard deviation of the distribution; (0, 1) is a random number distributed according to a standard normal distribution.

[0051] The acceptance probability is updated by calculating the distribution parameter vector using the likelihood ratio, as follows:

[0052] In the formula: L (·) is the sample likelihood function, and its mathematical definition is consistent with steps 1.2.2 and 1.2.3.

[0053] This method employs a uniform, uninformative prior and a symmetric normal proposal distribution. The ratio of the prior term to the proposal distribution term is always equal to 1; therefore, the acceptance probability only retains the likelihood ratio term. During calculation, a set of uniformly random numbers is generated in the interval [0, 1]. If the random number is less than... α acc Then the updated parameter is i prop Conversely, retain i curr .

[0054] The complete calculation process is demonstrated using a single iteration of 3-2 layer cohesion: setting the current distribution parameter vector. i curr =[2.825,0.210], perturbation generates candidate distribution parameter vectors i prop =[2.841,0.206], substituting the two sets of parameters into the likelihood function, we obtain L( i curr ) = 2.36 × 10 -12 L( i prop ) = 2.81 × 10 -12 ,calculate α acc =L( i prop ) / L( i curr Given 1.19, we take min(1,1.19)=1; generate a uniformly distributed random number of 0.623 in the interval 0~1. This value is less than the acceptance probability, so we accept the candidate parameters in this iteration and update... i curr = i prop .

[0055] The total number of global iterations is uniformly set to 50,000. The first 10,000 iterations are discarded as a pre-burn-in period, and the remaining valid posterior samples are counted. N mcmc =40,000 groups, and the average of the effective samples is taken as the posterior estimate of the parameters. These 40,000 posterior samples will be directly used for constructing the posterior prediction distribution and verifying the fitting effect in step 1.2.6, and the source of variance (standard deviation) of the covariance matrix in step 4.1. s The result is provided by the posterior estimation results of this step.

[0056] 1.2.6 Generate predicted samples and verify the goodness of fit Based on the sampling obtained from the MH-MCMC algorithm in step 1.2.5 N mcmc =40,000 valid posterior parameter samples were used to construct the posterior prediction distribution of key soil parameters. The posterior prediction distribution is in integral form over the posterior parameter samples:

[0057] In the formula: i t For the first t A group of valid posterior distribution parameter vector samples for the MH-MCMC algorithm; f opt This represents the optimal probability density distribution.

[0058] Each group in turn i t Substituting the optimal marginal distribution expression, batch-generated prediction samples of corresponding key soil parameters were obtained, and a complete posterior prediction sample set was obtained. By plotting histograms and kernel density curves, the distribution pattern of the predicted samples was compared with the original measured samples of layer 3-2 in Table 2. The sample distribution patterns highly overlapped, verifying that the optimal marginal distribution selected in this study has reliable fitting accuracy.

[0059] This step determines three key soil parameters: cohesion, internal friction angle, and compression modulus; it outputs the optimal marginal distribution type and posterior estimates (mean and standard deviation) of the distribution parameters for each parameter type, and transmits them to step S3 for probability integral transformation of the Copula model, and to step S4 for construction of the joint covariance matrix.

[0060] S2: Deterministic Finite Element Model Construction and Mesh Coordinate Extraction 2.1 Deterministic stratigraphic-structural finite element model construction The average values ​​of key soil parameters obtained from Table 2 were used as the fixed material parameters for soil layers to divide the geological regions. After configuring boundary constraints, a deterministic finite element numerical model was built. Based on the actual working conditions of the Wuhan Metro Line 12 shield tunnel project in this embodiment, a stratum-structure model was established. This model serves as the geometric carrier and mesh base for subsequent random field parameter assignment. The overall model diagram is shown below. Figure 3 As shown.

[0061] Based on the exploration report, the geometric information of the tunnel and strata was determined. The soil adopted the Mohr-Coulomb constitutive model, and the lining adopted the linear elastic constitutive model; corresponding displacement boundary constraints were set. Deterministic calculations were completed using the mean parameters of the soil layers, and the benchmark solution was output as a reference. The model contains 156,240 hexahedral solid elements. The number of elements determines the dimension of the covariance matrix in step S4 and controls the calculation scale.

[0062] 2.2 Batch extraction of unit and node information and processing of borehole observation data Using the deterministic numerical analysis model file built in this step, the script automatically identifies key fields of elements and nodes within the finite element input file, and batch exports all element numbers, element topology information (numbers of nodes associated with elements), and 3D coordinate information of all mesh nodes. The element topology information is the mesh mapping data showing the connection relationships between elements and nodes, used to record the node numbers associated with each element. The exported element topology information and 3D coordinate information of the mesh nodes are used to solve for the mesh center point coordinates in step 2.3, and also serve as the spatial index reference for constructing the random field covariance matrix in step S4.

[0063] The exploration borehole test data were compiled, coordinates and burial depths were matched, and measured samples with three-dimensional coordinates were generated as input for the joint ordinary kriging observation in step S4. The summary of the exploration measured data in this embodiment is shown in Table 5. These 16 sets of measured data serve as the joint ordinary kriging conditional correction observation constraints in step S4, ensuring that the parameter values ​​of the conditional three-dimensional conditional random field at the borehole locations are completely consistent with the field measured values.

[0064] Table 5 Summary of Exploration Data

[0065] 2.3 Calculation of grid center point coordinates Using the element topology information (element-node correspondence) and the three-dimensional coordinates of the mesh nodes exported in 2.2, the coordinates of the center point of each finite element mesh element are calculated; the element center point serves as the spatial positioning point for assigning random field parameters, and also as the spatial index reference for constructing the covariance matrix.

[0066] For conventional solid elements such as tetrahedral and hexahedral elements, the arithmetic mean of the coordinates of all nodes in the element is taken as the coordinates of the element center point:

[0067] In the formula: n n The number of nodes contained in the cell; x i , y i , z i ) is the unit number i The three-dimensional coordinates of each node, x c , y c , z c ) represents the three-dimensional coordinates of the center point of the finite element.

[0068] All N The coordinates of the 156,240 element center points are sorted according to the element number to generate a complete grid coordinate matrix, which is then fed into step S4 to complete the construction of the multi-parameter joint covariance matrix. After all geometric and coordinate data processing is completed, a basic Abaqus INP file is exported for later use; this INP file serves as the base file for assigning element-by-element values ​​to key soil parameters in the random field during batch automated modeling.

[0069] S3: Spatial autocorrelation function fitting and optimal Copula model fitting 3.1 Spatial Autocorrelation Fitting Based on the borehole geotechnical test and in-situ exploration data with accompanying three-dimensional spatial coordinates obtained in step S1, the horizontal autocorrelation distance is calculated using the autocorrelation function method. d x , d y Vertical autocorrelation distance d z .

[0070] 3.1.1 Preprocessing of raw borehole data Test data from all boreholes in the same soil layer were compiled and bound to three-dimensional spatial coordinates. For vertical correlation distance calculation, the data in the same layer were sorted from shallowest to deepest, and interpolated to generate a data sequence with equal depth intervals. The vertical sampling step size was uniformly set to Δz=0.25m, and the total number of measuring points in the depth sequence was determined. N seq Taken from the side expansion of the field flat shovel K D A total of 36 depth measurement points were measured in the experiment to eliminate the autocorrelation calculation error caused by the uneven spacing of the borehole sampling in the field.

[0071] 3.1.2 Detrending of Series and Calculation of Empirical Autocorrelation Function Polynomial fitting is used to remove the macroscopic trend term of parameter sequence changes with depth, resulting in a zero-mean stationary random sequence; the empirical autocorrelation function corresponding to each lag distance is calculated using the following formula:

[0072] In the formula: It is the empirical autocorrelation function; h j For the first j Spatial lag distance; N seq The total length of the standardized sequence; z ( i ) is the first i Depth measurement point parameters.

[0073] This function characterizes the spatial spacing. h j The strength of the linear correlation between two soil parameters reflects the law that the autocorrelation gradually decreases with increasing distance.

[0074] 3.1.3 Autocorrelation Function Model Fitting Three classical autocorrelation function models—squared exponential, single exponential, and cosine exponential—were selected. The Levenberg-Marquardt nonlinear least squares algorithm was used to fit the positive values ​​of the empirical autocorrelation function, fitting only the segments where the empirical autocorrelation was greater than 0, and a scaling parameter was established. b Autocorrelation distance d Conversion relationships: Squared exponential autocorrelation function (SQX):

[0075] Single exponential autocorrelation function (SNX):

[0076] Cosine exponential autocorrelation function (CSX):

[0077] In the formula: This is the theoretical autocorrelation function; h Spatial lag distance h j The abbreviation of .

[0078] 3.1.4 Evaluation of Fitting Performance and Determination of Optimal Correlation Distance Using the coefficient of determination R The formula for evaluating the fitting accuracy of the three types of autocorrelation models is as follows:

[0079] In the formula: n This represents the total number of hysteresis intervals involved in the fitting process. i ,k The lag distance is traversed and numbered, corresponding to the spatial lag distance. h j .

[0080] The autocorrelation distance corresponding to the maximum value of the coefficient of determination is selected as the optimal correlation distance; if the coefficients of determination of the three types of autocorrelation function models are all less than 0.1, it is determined that the spatial correlation of the soil layer parameters is weak, and the engineering experience values ​​are used as the correlation distances in the vertical and horizontal directions, respectively.

[0081] This embodiment uses the horizontal stress index of the lateral expansion of a flat spade. K D K generated from 36 depth measurement points in the experiment D The vertical correlation distance is calculated based on the depth sequence, and the curve of the horizontal stress index changing with depth is shown below. Figure 4 As shown in Table 4, the vertical correlation distance calculation results show that the single exponential autocorrelation function (SNX) has the best fitting accuracy (R²). 2 =0.93), corresponding to the vertical correlation distance. d z =2.94m. Using the exact same calculation process, the horizontal autocorrelation distance was calculated based on the borehole plane coordinates and in-situ horizontal test data, resulting in... d x =25.6m, d y =22.3m. The three sets of spatially correlated distances are supplied to step S4, which are directly used to construct the three-dimensional Gaussian covariance function and control the spatial decay rate of the random field.

[0082] Table 4 Calculation results of vertical correlation distance

[0083] 3.2 Fitting the Multi-parameter Copula Model Based on the optimal marginal distribution in step S1, the parameter samples are uniformly mapped in space and the sample interval [ε, 1-ε] is constrained. The maximum likelihood fitting of four types of Copula models is used, and the optimal model is selected by the AIC criterion and the correlation coefficient is output.

[0084] 3.2.1 Confirmation of Optimal Marginal Distribution The edge distributions of cohesion and internal friction angle were screened and verified using the AIC (Akaike Information Criterion).

[0085] Where: AIC K For the first k The numerical values ​​of the Akaike information criteria corresponding to the candidate distributions; For the first k One alternative distributionM k The probability density function; The complete set of distribution parameter vectors obtained from the maximum likelihood estimation; p k This represents the number of parameters to be estimated for the corresponding distribution model.

[0086] Select AIC K The distribution with the smallest value is taken as the optimal marginal distribution for fitting the Copula model. The optimal distribution for cohesion is log-normal, and the optimal distribution for internal friction angle is normal, which is consistent with the optimal distribution identified by the Bayesian model in step S1.2.

[0087] 3.2.2 Probability Integral Transform – Transformation to a Uniform Space Each set of key soil parameter samples is transformed by probability integral using the cumulative distribution function of its optimal edge distribution, mapping it to a uniform space to obtain the uniform edge samples required for Copula fitting. To avoid numerical singularities caused by boundary values ​​of 0 and 1, boundary truncation is performed on the transformation results. Calculation formula:

[0088] In the formula: F c (·) represents cohesion. c The corresponding optimal marginal cumulative distribution function; F φ (·) represents the internal friction angle. f The corresponding optimal marginal cumulative distribution function.

[0089] To avoid numerical instability caused by boundary values ​​of 0 or 1, the transformed values ​​are truncated at the boundary (taking...). e =1×10 -8 ),make sure u i , v i ∈[ε, 1-ε]. u i 、v i For the first i Uniform interval samples after mapping cohesion and internal friction angle; e The threshold is used for uniform spatial truncation.

[0090] 3.2.3 Substituting Multiple Types of Copula Functions into the Model Fitting For the 7 uniform samples after transformation in 3.2.2 ( u i , v iThe Gaussian, Frank, Gumbel, and t-Copula functions are fitted respectively. The joint distribution function expression of the four types of Copula is as follows:

[0091]

[0092]

[0093]

[0094] In the formula: u、v for u i 、v i Abbreviation; i For single-parameter Frank and Gumbel type Copula functions, the relevant scalar parameters control the degree of correlation between variables; r is the Pearson linear correlation coefficient, used in Gaussian and t-type Copula functions to describe the strength of the linear correlation between variables; n The degrees of freedom parameters of the Copula model corresponding to the t-class Copula function control the tail correlation of the distribution, and the degrees of freedom n The smaller the value, the stronger the tail correlation; Φ(·) is the standard normal distribution function, Φ -1 (·) is its inverse function (quantile function); t ν (·) represents the degrees of freedom. n of t Distribution function, t ν - ¹(·) is its inverse function; r , n , i Constructing the complete distribution parameter vector of the Copula model i Cop .

[0095] C (·) represents the Copula joint cumulative distribution function. c (·) represents the corresponding Copula probability density function, and the two satisfy the second-order partial derivative transformation relationship: .

[0096] Substituting the four types of Copula functions into their respective Copula models, we perform maximum likelihood estimation fitting, and the optimal parameter solution formula is as follows:

[0097] In the formula:c ( u i , v i ; i Cop ) represents the probability density function corresponding to the Copula function; i Cop This is the complete set of distributed parameter vectors for the Copula model (models corresponding to Gaussian, Frank, and Gumbel class Copula functions contain only single scalar parameters; models corresponding to t class Copula functions include correlation coefficients). r and degrees of freedom n (Two parameters).

[0098] After the optimal fit is completed, the log-likelihood value ln for each Copula model is calculated. L max It is used for model selection.

[0099] 3.2.4 Selection of the Optimal Copula Model The optimal Copula model is selected based on the log-likelihood value or the AIC (Akaike Information Criterion). The specific AIC calculation formula for the Copula model is as follows:

[0100] In the formula: ln L ( ) is the optimal parameter vector The maximum log-likelihood value obtained after fitting; This is the complete set of distribution parameter vectors for the optimal Copula model obtained by using the maximum likelihood estimation method; k This corresponds to the number of scalar parameters within the Copula model (the number of parameters in Gaussian, Frank, and Gumbel models). k =1, t-Copula includes the correlation coefficient. r Degrees of freedom n Two scalar parameters, k =2).

[0101] AIC Copula A smaller value indicates that the Copula model achieves a better balance between fitting accuracy and complexity. Selecting AIC... Copula The model corresponding to the minimum value is taken as the optimal Copula model. This example uses cohesion and internal friction angle as two sets of parameters, and compares the fitting results of four types of Copula models. Figure 5 As shown, the Copula model corresponding to the Gumbel class Copula function has the best fitting effect (the smallest AIC value), and the maximum log-likelihood ln L=13.80, corresponding to the Pearson linear correlation coefficient r =0.62, which can well characterize the upper tail correlation between cohesion and internal friction angle. The correlation coefficient of this optimal Copula model is output as the cross-covariance submatrix of the joint covariance matrix constructed in step 4.2. C cφ The amplitude coefficient directly determines the intensity of the cross-correlation between cohesion and internal friction angle.

[0102] 3.2.5 Generating random samples based on the optimal Copula model and inversely transforming them back to physical space A large number of uniform spatial sample pairs are generated from the fitted optimal Copula model, and then mapped back to the original physical parameter space using the optimal marginal distribution inverse cumulative distribution function (inverse CDF), generating soil parameter sample pairs that retain the original parameter cross-correlation structure; the sampling and inverse transformation formulas are as follows:

[0103]

[0104] In the formula: N sim The total number of simulated samples (in this embodiment, the value is fixed at 10,000 groups); The j-th group is a uniform interval random sample generated by the built-in Copula function of the optimal Copula model; This is the optimal Copula joint cumulative distribution function obtained through screening.

[0105] generated The sample pairs completely preserve the original c , f The cross-correlation structure between parameters can be directly used for subsequent random field modeling and stochastic finite element numerical analysis. The parameter samples generated by each Copula model are arranged according to cohesion. c (kPa), internal friction angle f (°) Save as a CSV file and output a model summary file synchronously, recording all Copula model fitting indices and optimal model parameter information.

[0106] S4: Construction of Joint Covariance Matrix and Generation of 3D Conditional Random Fields Combining the optimal edge distribution and parameter statistical characteristics obtained in step S1, and the autocorrelation distances in each direction calculated in step 3.1.4, d x =25.6m, d y =22.3m, d z=2.94m), step 3.2.4 obtains the cross-correlation coefficient of the parameters ( r =0.62), calculated in step 2.3 N =156240 finite element mesh center point coordinates, and 16 sets of three-dimensional measured parameter observation data obtained in step 2.2 (Table 5). Using the above inputs, a joint covariance matrix of multiple key soil parameters is constructed; using Ordinary Kriging interpolation combined with Karhunen-Loève (KL) expansion, a conditional three-dimensional conditional random field of key soil parameters is generated that simultaneously satisfies spatial autocorrelation, multi-parameter cross-correlation, and borehole measured data constraints.

[0107] 4.1 Define the covariance function The spatial autocorrelation structure of a single key soil parameter is characterized by a three-dimensional Gaussian covariance function, where any two spatial points... p 1( x 1, y 1, z 1) with p 2( x 2, y 2, z 2) The formula for calculating the autocovariance between them is:

[0108] In the formula: Δ x =| x 1- x 2|、Δ y =| y 1- y 2|、Δ z =| z 1- z 2| represents the difference in coordinates between two points along each coordinate axis. s ² represents the population variance of the parameter, derived from the standard deviation in Table 2. s ; d x , d y , d z They are respectively x , y , z The spatial autocorrelation distances in the three directions are taken from the calculation results in step 3.1.

[0109] For the cross-covariance between cohesion and internal friction angle, the spatial attenuation term maintains the same form as the autocovariance, and is calculated using the average of the relevant distances in the corresponding directions of the two parameters. The magnitude of the cross-covariance is determined by the cross-correlation coefficient of the parameters and the standard deviations of the two sets of parameters.

[0110] In the formula: r =0.62 is taken from the optimal Copula model fitting result in step 3.2.4; s c =3.57kPa and s φ =1.62° is taken from the standard deviation of the corresponding parameter in Table 2; , , This represents the average value of the autocorrelation distances in the corresponding directions for the two sets of parameters.

[0111] 4.2 Constructing the Joint Covariance Matrix Considering both cohesion and internal friction angle, two key soil parameters, a 2-dimensional model is constructed. N ×2 N joint covariance matrix C joint ; N Representing the total number of center points of the finite element mesh, the matrix is ​​divided into blocks as follows:

[0112] In the formula: C cc The autocovariance submatrix of cohesion (dimension) N × N ); C φφ The autocovariance submatrix of the internal friction angle (dimension) N × N ); C cφ The cross-covariance submatrix of cohesion and internal friction angle (dimensions) N × N ).

[0113] Joint covariance matrix C joint It fully characterizes the spatial autocorrelation and cross-correlation structure of the two types of parameters at all grid points.

[0114] 4.3 Constructing the covariance matrix between observation points Based on the three-dimensional measured key soil parameter observation samples compiled in step 2.2 (Table 5 contains a total of 16 sets of borehole observation data), a system was constructed. n obs × n obs observation point covariance matrix C oo Each element in the matrix is ​​calculated using the autocovariance function and crosscovariance function defined in step 4.1, according to its parameter category. To improve the numerical stability of matrix inversion, a small perturbation term of 1×10⁻⁶ is superimposed on all elements along the main diagonal of the matrix.-10 .

[0115] Simultaneously construct 2 N × n obs Grid point-observation point cross-covariance matrix C uo It is used to characterize the covariance relationship between the center points of all finite element mesh elements and the observation points of each measured borehole.

[0116] 4.4 Combined Ordinary Kriging Interpolation – Conditional Mean Estimation The conditional average of a three-dimensional conditional random field is calculated using the joint ordinary kriging method. The matrix form of the joint ordinary kriging equations is as follows:

[0117] The weight matrix is ​​calculated. w (dimension) n obs ×2 N Then, the conditional average vector is calculated:

[0118] In the formula: y obs This is the vector of measured observations from Table 5 in step 2.2. m Krig It is a dimension of 2 N A 1×1 average vector containing the joint ordinary kriging estimate of cohesion and internal friction angle at all grid points.

[0119] Conditional average m Krig This represents the optimal linear unbiased estimation result of the spatial distribution of key soil parameters under the constraint of known borehole observation data.

[0120] 4.5 Karhunen-Loève Expanding – Unconditional Sample Generation For the joint covariance matrix C joint Perform eigenvalue decomposition:

[0121] In the formula: V The eigenvector matrix; Λ=diag( l 1, l 2, ..., l 2N ) is an eigenvalue diagonal matrix.

[0122] The eigenvalues ​​are sorted from largest to smallest. The eigenvalues ​​are truncated according to the energy ratio criterion of a cumulative energy percentage of at least 0.95, retaining the principal eigenmode. Karhunen-Loève expansion is used to generate unconditional random samples. Discriminant formula:

[0123] In this embodiment, the energy percentage threshold is fixed at 0.95, which significantly reduces the random dimension while ensuring the fitting accuracy of unconditional random samples (in this embodiment). N =156240, after dimensionality reduction n KL (Approximately 1843 orders, a reduction of 99.4%). This truncation process ensures that the unconditional random samples generated in step 4.6 retain more than 95% of the spatial variation information, while significantly reducing the storage and computation required for subsequent batch INP generation.

[0124] Unconditional random samples are calculated based on the truncated feature modes, using the following formula:

[0125] In the formula: g ~ (0, I )for n KL A standard normal random vector.

[0126] This method, known as Karhunen-Loève (KL) expansion, can efficiently generate a large number of unconditional random samples while fully preserving the spatial structure of the covariance of soil and rock parameters.

[0127] 4.6 Conditional Modification – Implementation of Conditional Random Fields By combining ordinary kriging conditional correction, unconditional random samples are transformed into three-dimensional conditional random field samples, ensuring that the samples are accurately consistent with the measured data at the observation points. The complete correction steps are as follows: ① Match the nearest neighbor index position of each measured observation point on the grid; ② Extract the values ​​of unconditional random samples at the observation points. ψ obs ; ③ Calculate the joint ordinary kriging estimate of the unconditional sample across all grid points using the joint ordinary kriging weights. ψ Krig :

[0128] ④ The following formula is used to complete the condition correction:

[0129] The corrected 3D conditional random field sample perfectly matches the measured parameters at the borehole observation points, while fully preserving the original spatial autocorrelation and multi-parameter cross-correlation structures of the parameters. The 3D distribution of the soil and rock parameters generated in this embodiment is as follows: Figure 6 As shown. By repeating the entire sampling, interpolation, and condition correction process, multiple sets of independent three-dimensional conditional random field implementation samples can be generated for subsequent stochastic finite element numerical analysis.

[0130] Thus far, all research results from S1 to S4 include: the sample statistics in Table 2 of step 1.1, the optimal marginal distribution in step 1.2, the spatial correlation distance in step 3.1.4 (Table 4), and the Copula correlation coefficient in step 3.2.4. r =0.62, 200 sets of three-dimensional conditional random fields in step 4.6. Based on the deterministic numerical analysis model (N=156240) built in step 2.1, and combined with the 200 sets of conditional key soil random parameters generated in step 4.6, the uncertainty analysis of geotechnical parameters is carried out using the non-intrusive stochastic finite element method. The entire calculation process is executed in a closed loop by a Python automated script, without the need for manual editing of the INP file.

[0131] The script automatically reads the element-parameter correspondence data table output in step 4.6 and batch modifies the INP input files: it iterates through all N=156240 soil solid elements in the model, creating an independent element set for each mesh and assigning a dedicated material entry; it reads key soil random parameters at the element scale and assigns cohesion, internal friction angle, and compression modulus to the corresponding element material properties. After assigning values ​​to a single set of 3D conditional random fields, the file is automatically saved. By iteratively generating 200 sets of INP calculation files with independent parameters, the script eliminates repetitive manual work and human error.

[0132] The 200 batch-generated INP files were submitted to Abaqus for parallel numerical computation to calculate engineering indicators such as ground settlement and lining internal forces in the shield tunnel. The average value, standard deviation, and 95% two-sided confidence intervals of the maximum surface settlement and maximum lining internal forces were extracted and statistically analyzed. The statistical results can quantitatively characterize the impact of the spatial variability of key soil parameters on the tunnel structure response, forming a complete technical chain of "dispersion degree of key soil parameters - spatial variability characteristics of random fields - fluctuation range of tunnel structure response," providing quantitative support for tunnel structure optimization design and construction risk prediction.

[0133] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for constructing and analyzing random fields considering uncertainties in geotechnical parameters, characterized in that, Includes the following steps: S1: Collect three-dimensional coordinate data and physical and mechanical parameter data of each soil layer in the site, divide the soil layers according to depth and collect them in layers, calculate the statistical characteristics of the parameters of each layer, determine the key soil parameters, and determine the optimal marginal distribution of each key soil parameter and the posterior estimate of its distribution parameters based on Bayesian inference. S2: Construct a deterministic finite element numerical model using the average values ​​of physical and mechanical parameters, extract the element topology information and three-dimensional coordinate information of the mesh nodes of the model, and calculate the three-dimensional center point coordinates of each mesh element; S3: Construct a depth data sequence using the three-dimensional coordinate data from step S1, calculate the empirical autocorrelation value, fit the autocorrelation function to determine the spatial correlation distance in each direction, and combine the optimal edge distribution determined in step S1 to fit the Copula model to obtain the correlation coefficient between parameters. S4: Based on the coordinates of the three-dimensional center point of the grid cell, the spatial correlation distance, the correlation coefficient of the optimal Copula model, and the distribution parameters estimated from the posterior, construct the joint covariance matrix of the key soil parameters, use Karhunen-Loève expansion to generate unconditional random samples, combine the physical and mechanical parameter data from step S1 to perform condition correction, and output the three-dimensional conditional random field of the key soil parameters.

2. The method according to claim 1, characterized in that, The physical and mechanical parameter data in step S1 are obtained through drilling sampling, indoor geotechnical tests, in-situ detection, and wave velocity testing. The in-situ detection and wave velocity test data both have continuous depth coordinates. The statistical characteristics include maximum value, minimum value, average value, standard deviation, and coefficient of variation. The key soil parameters are those with larger coefficients of variation among the physical and mechanical parameters of each soil layer.

3. The method according to claim 1, characterized in that, The specific process of Bayesian inference and posterior estimation in step S1 is as follows: Four candidate probability distribution models are selected, and a uniform prior distribution is set for the distribution parameters of each model; the likelihood function of each model is constructed based on the data of key soil parameters; according to Bayes' theorem, the posterior distribution of the distribution parameters is calculated from the prior distribution and the likelihood function; the model evidence of each model is calculated by Monte Carlo integration, and the posterior probability of each model is calculated based on the model evidence, and the model with the largest posterior probability is selected as the optimal marginal distribution; the distribution parameters of the optimal marginal distribution are posteriorly sampled using the Metropolis-Hastings Markov chain Monte Carlo method to obtain the posterior estimate of the distribution parameters, and a prediction sample is generated by the posterior estimate, and the fit of the optimal marginal distribution is verified by comparing it with the measured data.

4. The method according to claim 3, characterized in that, The alternative probability distribution models are normal distribution, log-normal distribution, gamma distribution and Weibull distribution. All four types of distributions include location parameters and scale parameters, while gamma and Weibull distributions additionally include shape parameters.

5. The method according to claim 4, characterized in that, The uniform prior interval of the position parameter is [0, 2]. μ The uniform prior interval for scale and shape parameters is [0.1]. σ ,2 σ ], μ This is the average value. σ To eliminate the standard deviation, when calculating model evidence, a logarithmic transformation is used to subtract the maximum value to eliminate the underflow of Monte Carlo integral values.

6. The method according to claim 1, characterized in that, In step S2, when constructing the deterministic finite element numerical model, the average value of the physical and mechanical parameters is used as the fixed parameter, and boundary conditions and loads are set. The element topology information and the three-dimensional coordinate information of the grid nodes of the model are extracted. The element topology information is the grid mapping data of the connection relationship between elements and nodes. The associated nodes corresponding to each grid element are retrieved using the element topology information, and the three-dimensional coordinate information of the grid nodes corresponding to the associated nodes is read and the average value is calculated. The three-dimensional center point coordinates of the grid element are thus calculated.

7. The method according to claim 1, characterized in that, In step S3, the specific process of autocorrelation function fitting and Copula model fitting is as follows: determine the equal-interval sampling step size, calculate the empirical autocorrelation function values ​​of each physical and mechanical parameter under different spatial lag distances; use multiple autocorrelation function models for fitting, select the optimal autocorrelation function model based on the goodness of fit, calculate the spatial correlation distance of each soil layer's physical and mechanical parameters in the vertical and horizontal directions; use multiple Copula models to fit the correlation structure between each key soil parameter, select the optimal Copula model based on the fitting index, and output its correlation coefficient.

8. The method according to claim 7, characterized in that, The autocorrelation function model selected includes three types: quadratic exponential, single exponential, and cosine exponential. The positive values ​​of the empirical autocorrelation function are fitted using the Levenberg-Marquardt nonlinear least squares algorithm. The goodness of fit of the three models is considered when... R 2 When the value is less than 0.1, the relevant distances in the vertical and horizontal directions are determined by engineering experience.

9. The method according to claim 7, characterized in that, The process of fitting the optimal Copula model is as follows: Based on the optimal edge distribution and its distribution parameters determined in step S1, the cumulative distribution function corresponding to the optimal edge distribution is obtained, and the data of each key soil parameter are mapped to a uniform space through the cumulative distribution function; the mapped data are then adjusted to the boundary so that they are located in […]. ε ,1- ε Within the interval, ε =1×10 -8 ; Four types of Copula models—Gaussian, Frank, Gumbel, and t-Copula—were used to fit the correlation structure between key soil parameters using maximum likelihood estimation. The AIC (Akaike Information Criterion) was used to validate the models, and the Copula model with the smallest AIC value was selected as the optimal Copula model. Its correlation coefficient was then output.

10. The method according to claim 1, characterized in that, Step S4 specifically involves: constructing an autocovariance matrix based on the coordinates of the grid cell center points and the spatial correlation distances in each direction; determining the amplitude using the posterior-estimated distribution parameters; constructing the cross-covariance matrix between parameters using the optimal Copula correlation coefficient; and forming the joint covariance matrix of the key soil parameters. Its eigenvalues ​​are decomposed, and the eigenvalues ​​are truncated to a cumulative characteristic energy ratio of not less than 0.95 while retaining the principal characteristic modes. Karhunen-Loève expansion is used to generate unconditional random samples. Combined with joint ordinary kriging and physical and mechanical parameter data, conditional corrections are performed to output a three-dimensional conditional random field of key soil parameters.

Images