CatBoost–Bayesian Hybrid Model Adaptively Coupled with Modified Theoretical Equations for Estimating the Undrained Shear Strength of Clay

Abstract

The undrained shear strength of clay is an important index for the calculation of the bearing capacity of the foundation soil, the calculation of the soil pressure of the foundation pit, and the analysis of the slope stability. Therefore, the purpose of this paper is to conduct a comprehensive study of the combined use of machine learning with clay theoretical equations to estimate it. Under the Bayesian framework, the CatBoost algorithm (CatBoost–Bayesian) based on Bayesian optimization algorithm was developed to obtain the feature importance level of soil parameters affecting the undrained shear strength of clay, so as to adaptively couple the theoretical equation of undrained shear strength of $K_0$ consolidated clay, which was derived from the modified Cambridge model. Then, the theoretical equation of undrained shear strength of the isotropically consolidated clay was established from the critical state of the clay parameters. Finally, it was illustrated and verified using the experimental samples of Finnish clay. The results indicate that the theoretical equation established by the overconsolidation ratio and effective overburden pressure parameters can well estimate the undrained shear strength of isotropically consolidated clays, and the parameter uncertainty can be considered explicitly and rigorously.

1. Introduction

The composition and physical and mechanical properties of clayey soil are important indexes in geotechnical engineering, among which the undrained shear strength index of clay is the most important index for calculation of the bearing capacity of foundation soil, calculation of the soil pressure of the foundation pit, and analysis of slope stability. Cross plate shear tests in the field, direct shear tests in the laboratory, triaxial compression tests, and other such direct measurement methods are time-consuming and expensive [1]. To evaluate through indirect methods, theoretical derivation and empirical methods are effective methods. Based on many experiments, Mesri and Ladd et al., respectively, proposed empirical equations for the undrained shear strength of clay [2,3]. Jiang Shuihua et al. simulated the uncertainty of soil parameters and proposed an effective non-stationary random field model for undrained shear strength parameters [4]. Indirect methods also include measurements obtained from dilatometer tests (DMT) based on empirical equations [5,6].

The application of machine learning in the field of geotechnical engineering has gradually become a research hotspot, including slope stability testing, TBM performance evaluation, rockburst vibration estimation, and pile foundation evaluation, etc., and has been proposed, expanded upon and applied by many researchers [7,8,9,10]. The emergence of artificial intelligence techniques holds great potential for solving soil parameter estimation problems involving complex soil–structure interactions [11,12,13]. However, compared with other directions in the field of geotechnical engineering, research devoted to soil parameter estimation is still very limited. The CatBoost algorithm can effectively solve gradient offset and prediction offset problems, improving the accuracy and generalization ability of the calculation [14]. Therefore, it has been widely used in crop evapotranspiration estimation in hydrology, data-driven seismic performance evaluation of corroded RC columns in civil engineering, and groundwater salinity prediction in multi-layer coastal aquifers in geotechnical engineering [14,15,16,17,18].

Moreover, the combination of Bayesian optimization with machine learning algorithms is still relatively rare, and Zhang et al. use extreme gradient boosting and random forest based on Bayesian optimization in geotechnical engineering to effectively capture the relationship between undrained shear strength and various fundamental soil parameters [19]. Ho et al. used the good performance of Bayesian, functional, and meta-ensemble machine learning models to generate land subsidence susceptibility (LSS) maps [20]. The Cambridge model and modified Cambridge model are two elastic–plastic constitutive models commonly used in geotechnical engineering to describe the shear properties of soils [21]. They are typical of the constitutive relationships of soils under undrained conditions. The Cambridge model describes the shear properties of soils under circular stress paths, while the modified Cambridge model modifies the stress paths in the plastic potential function to better describe the shear properties of soils under non-circular stress paths [22].

Under the Bayesian framework, the CatBoost algorithm (CatBoost–Bayesian) based on a Bayesian optimization algorithm was developed to obtain the feature importance level of soil parameters affecting the undrained shear strength of clay by using the experimental samples of Finnish clay, so as to adaptively couple the theoretical equation of undrained shear strength of $K_0$ consolidated clay, which was derived from the modified Cambridge model. Then, the theoretical equation of undrained shear strength of isotropically consolidated clay was established from the critical state of clay parameters. Finally, this paper analyzed the uncertainty of the parameters of the equation and verified the rationality of the calculation results. The rest of this study was arranged as follows. In Section 2, the CatBoost algorithm, the principle of Bayesian optimization algorithm based on random forest and K-fold cross-validation, as well as the theoretical equation derived from the revised Cambridge model, and the selection of quantitative evaluation indicators are introduced, respectively. The Finnish clay database used for training and validating the CatBoost–Bayesian hybrid model and theoretical equations is presented in Section 3, and a sensitivity analysis of different types of clay was performed to show the results of the importance of clay parameters. Finally, various verifications of the calculation results of the undrained shear strength theoretical equation of isotropically consolidated clay were carried out, and the results of the similar mixed model of CatBoost–Bayesian were compared and analyzed. Section 4 discusses the conclusions of this study, followed by references.

2. Materials and Methods

2.1. CatBoost Algorithm

The CatBoost algorithm was developed by researchers and engineers at Yandex, the largest Russian search engine company, in 2017, and opened source code in April [14,23]. It is the latest algorithm in the open-source field in the world today, and surpasses the XGBoost algorithm and the LightGBM algorithm in terms of performance. The name comes from “Category” and “Boosting”, and belongs to the family of Boosting algorithms. The CatBoost algorithm is an improvement in the framework of the GBDT algorithm. It effectively solves the problem of gradient bias and prediction shift, avoids the occurrence of overfitting, and improves calculation accuracy and generalization ability (Figure 1). The details are as follows.

• The CatBoost algorithm can handle categorical features in GBDT features better, and the simplest way is to use the average value of the corresponding labels to replace them. In the decision tree, the label average value will be used as the criterion for node splitting. This method is known as greedy target-based statistics, or greedy TS. However, this method has obvious drawbacks, so greedy TS is improved by adding prior distribution terms to reduce the effect of noise and low frequency categorical data on the data distribution [14,23].

$${\widehat{x}}_k^i=\frac{{{\displaystyle \sum }}_{j=1}^{p-1}\left[x_{{\sigma }_{j,k}}=x_{{\sigma }_{p,k}}\right]\cdot Y_{{\sigma }_j}+a\cdot p}{{{\displaystyle \sum }}_{j=1}^{p-1}\left[x_{{\sigma }_{j,k}}=x_{{\sigma }_{p,k}}\right]+a}$$

(1)

where p is the added prior term and a is usually a weighting factor greater than 0. For regression problems, the prior term can be taken as the average value of the dataset in general.

2.

The prediction shift is caused by the gradient bias. To overcome this problem, CatBoost proposes a new algorithm called ordered boosting (Algorithm 1).

Algorithm 1: Ordered boosting pseudo-code algorithm

input: ${\left\{\left(X_k,Y_k\right)\right\}}_{k=1}^n,$ $I;$

$\sigma \leftarrow$random permutation of $\left[1,n\right];$

$M_i\leftarrow 0$ for $i=1\dots n;$

for $t\leftarrow 1$ to $I$ do

for $i\leftarrow 1$ to $n$ do

$r_i\leftarrow y_i-M_{\sigma \left(i\right)-1}\left(X_i\right);$

for $i\leftarrow 1$ to $n$ do

$\mathsf{\Delta }M\leftarrow$learn model $\left(\left(X_j,r_j\right):\sigma \left(j\right)\le i\right);$

$M_i\leftarrow M_i+\mathsf{\Delta }M;$

From the Algorithm 1, $X_i$ denotes a sample, $M_i$ denotes a separate model, and model $M_i$ is obtained by training with a training set that does not contain sample $X_i$.

2.2. Bayesian Optimization Algorithm (SMAC) and k-Fold Cross-Validation

Hyperparameter optimization is a combinatorial optimization problem that cannot be optimized by gradient descent as general parameters. Evaluating a set of hyperparameter problems often involves issues such as resource allocation and efficiency. The simpler hyperparameter optimization methods include manual search, grid search, and random search. In this study, we used the currently popular adaptive hyperparameter search method of Bayesian optimization [24,25,26]. Based on the combinations of hyperparameters that were already tested, this method can predict the next combination that is more likely to bring the greatest benefit. This study is implemented using the optuna automatic hyperparameter optimization framework.

SMBO stands for sequential model-based optimization. The so-called serialization refers to the optimization by iterative method one trial at a time. SMBO is a specific implementation form of Bayesian optimization.

SMAC, proposed by Hutter et al. [27], stands for sequential model-based optimization for general algorithm configuration; the model originates from the random forest model. The random forest (RF) algorithm is a combination of the Bootstrap Aggregating algorithm proposed by Breiman in 1996 and the stochastic subspace algorithm proposed by Ho [28] in 1998. SMAC was initially designed to remove some of the limitations of SMBO to make it applicable to general algorithm configuration problems with many classification parameters and benchmark sets of instances, and to solve the case in which the parameter type cannot be discrete in the Gaussian regression process.

During the process of model training, the problem of data overfitting often occurs. The model can match the training data well, but cannot predict the data outside the training set well. This affects the accuracy of the final evaluation results.

K-fold cross-validation (CV) is widely used as a criterion for model selection [29], and was originally created to reduce the computational cost of leave-one-out cross-validation (LOOCV) [30]. It divides a dataset into K equal samples, of which K-1 samples are used to build the model for training, and the remaining 1 sample is used for validation. During this K-time iteration, each part of K is successively assigned as validation data.

Jung [31] proposed a new method of choosing K-fold to ensure K = ln(n) and guarantee $n/\mathrm{K} > 3d$ simultaneously. After verifying the feasibility of it, where n denotes a data point and d denotes the number of parameters (in this study, n = 202, d = 12, K = 5.3), and combining it with Pham [32] method, K was setted to 5. Therefore, this study used a five-fold cross-validation method for validation (Figure 2).

To reduce the reliance on empirical rules and inefficient brute force search, Bayesian optimization methods are applied to determine the appropriate model hyperparameters for the CatBoost algorithm (CatBoost–Bayesian). With appropriate model hyperparameters, the CatBoost–Bayesian hybrid model is able to provide a data-driven ranking of clay feature importance and properly uncover the intrinsic linkage of soil parameter essences; however, it cannot give full interpretability. To address this issue, the clay feature importance is adaptively coupled with the theoretical equations derived from the modified Cambridge model to estimate the undrained shear strength of isotropically consolidated clay under a Bayesian probabilistic framework, taking into account the parameter and model uncertainties, which are explained in detail in Section 2.3 and Section 2.4.

2.3. Theoretical Equation

Deriving the undrained shear strength ratio $S_u/{\sigma }_{vi}^{’}$ from the soil intrinsic relationship is a common method for solving soil strength. Based on the Cambridge model, Ohta et al. [33] proposed the undrained shear strength derivation equation under the condition of triaxial $K_0$ ($K_0$ is the coefficient of earth pressure at rest) consolidation.

$${\left(\frac{S_u}{{\sigma }_{vi}^{\prime }}\right)}_{\mathrm{OCA}}=OC{R}^{\stackrel{-}{\Lambda }}\cdot {\left(\frac{S_u}{{\sigma }_{v0}^{\prime }}\right)}_{\mathrm{NCA}}$$

(2)

where OCA is the overconsolidated state of anisotropic clay, NCA is the normally consolidated state of anisotropic clay, the overconsolidation ratio is $OCR={\sigma }_{v0}^{\prime }/{\sigma }_{vi}^{\prime }$, $\stackrel{-}{\Lambda }=1-\stackrel{-}{\kappa }/\lambda$, $\stackrel{-}{\kappa }$ is the slope of the swelling lines described by the effective overburden pressure, and $\lambda$ is the slope of the anisotropical normally consolidated soil compression curve:

$$\stackrel{-}{\kappa }=\frac{\left(e_i-e_0\right)}{\left(ln{\sigma }_{v0}^{\prime }-ln{\sigma }_{vi}^{\prime }\right)}$$

(3)

$$\kappa =\frac{\left(e_i-e_0\right)}{\left(ln{p}_0^{\prime }-ln{p}_i^{\prime }\right)}$$

(4)

where $\kappa$ is the slope of the swelling lines described by the effective average stress, $p_0^{\prime }$ is the effective average vertical pressure in the anisotropical normally consolidated state, and $p_i^{\prime }$ is the effective average overburden pressure. $e_i$ is the void ratio, and $e_0$ is the initial void ratio.

$$\frac{p_0^{\prime }}{p_i^{\prime }}=OC{R}^{\left(\stackrel{-}{\kappa }/\kappa \right)}$$

(5)

Substitute $p_0^{\prime }=\frac{1}{3}\left(1+2K_{0nc}\right)\cdot {\sigma }_{v0}^{\prime }$ and $p_i^{\prime }=\frac{1}{3}\left(1+2K_0\right)\cdot {\sigma }_{vi}^{\prime }$ into Equation (5), to obtain the following:

$$\stackrel{-}{\kappa }=\kappa \cdot \frac{ln\left(\frac{1+2K_{0nc}}{1+2K_0}\right)+ln\left(OCR\right)}{ln\left(OCR\right)}$$

(6)

where $K_0$ is the coefficient of earth pressure at rest in the soil in the overconsolidated state, $K_{0nc}$ is the coefficient of earth pressure at rest in the soil in the normally consolidated state, and the undrained shear strength ratio of the normally consolidated soil is

$${\left(\frac{S_u}{{\sigma }_{v0}^{\prime }}\right)}_{\mathrm{NCA}}=\frac{1+2K_0}{3}\frac{M}{2}exp\left(-\Lambda \pm \frac{\Lambda }{M}{\eta }_0\right)$$

(7)

where $\Lambda =1-\kappa /\lambda$, ${\eta }_0=3\cdot \left(1-K_0\right)/\left(1+2\cdot K_0\right)$, $M$ is the critical state stress ratio, $M=\mathrm{q}/p$, $q$ is generalized shear stress, $q=\left({\sigma }_1-{\sigma }_3\right)$, $p$ is the average principal stress, and $p=\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)/3$. When ${\eta }_0=0$ and $K_0=1$, the ratio of undrained shear strength of the isotropical normally consolidated state and isotropically overconsolidated state clay is obtained, respectively.

$${\left(\frac{S_u}{{\sigma }_{v0}^{\prime }}\right)}_{\mathrm{NCI}}=\frac{M}{2^{\Lambda +1}}$$

(8)

$${\left(\frac{S_u}{{\sigma }_{vi}^{\prime }}\right)}_{\mathrm{OCI}}=OC{R}^{\Lambda }\cdot {\left(\frac{S_u}{{\sigma }_{v0}^{\prime }}\right)}_{\mathrm{NCI}}=OC{R}^{\Lambda }\frac{M}{2^{\Lambda +1}}$$

(9)

where OCI is the overconsolidated state of isotropic clay, and NCI is the normally consolidated state of isotropic clay. Based on the experimental results of isotropic consolidation and swelling tests, Karube [34] proposed an empirical equation $M=1.75\cdot \Lambda$.

From Equations (8) and (9), the equations for calculating the undrained shear strength of the isotropical normally consolidated and isotropically overconsolidated clays can be obtained.

$${\left(S_u\right)}_{\mathrm{NCI}}={\sigma }_{v0}^{\prime }\cdot \frac{1.75\cdot \Lambda }{2^{\Lambda +1}}={\sigma }_{v0}^{\prime }\cdot \frac{1.75\cdot \left(1-\frac{\kappa }{\lambda }\right)}{2^{\left(2-\frac{\kappa }{\lambda }\right)}}$$

(10)

$${\left(S_u\right)}_{\mathrm{OCI}}={\sigma }_{vi}^{\prime }\cdot OC{R}^{\Lambda }\cdot \frac{1.75\cdot \Lambda }{2^{\Lambda +1}}={\sigma }_{vi}^{\prime }\cdot OC{R}^{\left(1-\frac{\kappa }{\lambda }\right)}\cdot \frac{1.75\cdot \left(1-\frac{\kappa }{\lambda }\right)}{2^{\left(2-\frac{\kappa }{\lambda }\right)}}$$

(11)

In the Cambridge model, in isotropically consolidated clays, $\lambda =0.434\cdot Cc$, $\kappa =0.434\cdot Cs$ [35].The initial void ratio can be solved by equation $e_0=2.5\cdot C_c+0.25$ [36]. Ladd [37] proposed an empirical equation for computing the undrained shear strength of clay based on the results of indoor and outdoor tests and examples of foundation failure.

$$\frac{S_u}{{\sigma }_{vi}^{\prime }}=S\cdot {\left(OCR\right)}^m$$

(12)

where $S$ is under normal consolidation conditions, $m$ is the strength growth index, and $m$ is affected by $1-\kappa /\lambda$ or $1-C_s/C_c$ according to the critical state theory [21].

2.4. Bayesian Perspective of Unified Undrained Shear Strength Equation

Unlike classical statistics, which is based on the frequentist approach, the CatBoost–Bayesian hybrid model is based on Bayesian statistics and machine learning. The CatBoost–Bayesian hybrid model uses Bayesian theory to combine engineering judgment and empirical prior knowledge with data to derive posterior distributions. In geotechnical analysis and prediction, engineers usually consider data from multiple sources [38]. This study updates the intrinsic model parameters (posterior distribution) with information from laboratory experimental data, field test data, and engineering experience [39]. Considering the differences in soil area and experimental methods, according to Equations (7)–(12), the effective overlying pressure ${\mathsf{\sigma }}_{{\mathrm{v}}_{\mathrm{i}}}^{\prime }$ is used together, $\mathrm{and} {\mathsf{\sigma }}_{\mathrm{vi}}^{\prime }$ and ${\mathsf{\sigma }}_{\mathrm{v}0 }^{\prime }$ are converted by the correction coefficient $\mathsf{\beta }$ to complete the unification of the undrained shear strength equations of isotropic normal consolidation and isotropic overconsolidation, derived from the constitutive model.

$$S_u=f\left(\xi \right)+\epsilon$$

(13)

where $\xi =\left[\alpha ,\beta ,m,{\sigma }_{v_i}^{\prime },OCR\right]$ represents the constitutive model parameter, $f\left(\xi \right)=\alpha {\left({\sigma }_{v_i}^{\prime }\right)}^{\beta }{\left(OCR\right)}^m$, $\alpha =1.75m/2^{\left(1+m\right)}$, $n=1-k/\lambda =1-c_s/c_c$, and $\epsilon$ is a Gaussian random variable with mean $u_{\epsilon }$=0 and standard deviation ${\sigma }_{\epsilon }$= 0.215 [40]. Conditional probability theory is used to update the $S_u$ posterior probability density function (PDF) based on prior information and laboratory test data of random variables. The calculation is as follows:

$$p\left(\xi |S_u\right)=KL\left(S_u|\xi \right)p\left(\xi \right)$$

(14)

where K is the normalization constant; $p\left(\xi \right)$ is the prior distribution of the key constitutive parameter $\xi$, and $L\left(S_u|\xi \right)$ is the likelihood function reflecting the probability relationship between $\mathrm{the} S_u$ of the laboratory test (site test) and constitutive model parameters.

2.4.1. Prior Distribution

Information about model parameters is usually limited given laboratory/field tests. Therefore, non-informative prior knowledge is used in this study. Assuming that a single random variable is independent, the prior distribution is equal to the product of the prior distributions of all constitutive model parameters. In order that there is no preference for any value in the possible range of parameters, uniform prior distribution is usually used. Uniform prior distribution is used to reflect the knowledge state of the model parameters in the Bayesian framework [41], expressed as

$$p\left(\xi \right)={\displaystyle \prod }_{i=1}^{n}p\left({\xi }_i\right)$$

(15)

$$p\left(\xi \right)=\{\begin{array}{cc}\frac{1}{{\mu }_{imax}-{\mu }_{imin}}·\frac{1}{{\sigma }_{imax}-{\sigma }_{imin}}& {\mu }_iϵ\left[{\mu }_{imin},{\mu }_{imax}\right],{\sigma }_iϵ\left[{\sigma }_{imin},{\sigma }_{imax}\right]\\ 0& others\end{array}.$$

(16)

where $n$ is the number of constitutive model parameters, and ${\mu }_{imin}$, ${\mu }_{imax}$, ${\sigma }_{imin}$ and ${\sigma }_{imax}$ are the lower and upper limits of ${\mu }_i$ and ${\sigma }_i$ respectively.

2.4.2. Likelihood Function

The likelihood function reflects the degree of consistency between constitutive model parameter A and the statistical characteristics of laboratory test/field test b, which can be expressed by multivariate normal distribution function [42]:

$$L\left(S_u|\xi \right)={\left(\frac{1}{\sqrt{2\pi }{\sigma }_{S_u}}\right)}^N{\displaystyle \prod }_{i=1}^N\exp \left[\frac{{\left(f_i\left(\xi \right)-{\mu }_{S_u}\right)}^2}{2{\sigma }_{S_u}^2}\right]$$

(17)

where N is the number of laboratory tests/field tests. $\mu S_u$ and $\sigma S_u$ are the mean and standard deviation of $S_u$ respectively, and $f_i\left(\xi \right)$ is the value estimated by the parameters of the constitutive model. With the continuous development of the field of geotechnical engineering driven by data, it is less difficult to obtain prior information than in previous engineering practices. In order to obtain updated posterior $\xi$, Markov chain Monte Carlo simulation (MCMC) is used to obtain discrete samples based on prior distribution and likelihood function simulation [43,44]. The basic idea of MCMC simulation is to repeatedly extract samples from the target distribution by converging with the Markov chain of the target distribution.

2.4.3. Uncertainty Analysis of the Unified Constitutive Model

The uncertainty of the undrained shear strength equations of isotropic normal consolidation and isotropic overconsolidation after updating the constitutive model parameters is studied by using a first-order approximation method. The mean value and standard deviation of a can be approximated as [45,46]

$${\mu }_{S_u}=\alpha {\left({\mu }_{{\sigma }_{vi }^{\text{'}}}\right)}^{\beta }{\left({\mu }_{OCR}\right)}^m$$

(18)

$$\begin{array}{c}{\sigma }_{S_u}^2={\left(\frac{\partial S_u}{\partial {\sigma }_{vi }^{\prime }}\right)}^2{\sigma }_{{\sigma }_{vi }^{\prime }}^2+{\left(\frac{\partial S_u}{\partial OCR}\right)}^2{\sigma }_{OCR}^2+{\sigma }_{\epsilon }^2\\ \approx {\left[\frac{S_u\left({\mu }_{{\sigma }_{vi }^{\prime }}+{\sigma }_{{\sigma }_{vi }^{\prime }}\right)-S_u\left({\mu }_{{\sigma }_{vi }^{\prime }}-{\sigma }_{{\sigma }_{vi }^{\prime }}\right)}{2}\right]}^2+{\left[\frac{S_u\left({\mu }_{OCR}+{\sigma }_{OCR}\right)-S_u\left({\mu }_{OCR}-{\sigma }_{OCR}\right)}{2}\right]}^2+{\sigma }_{\epsilon }^2\end{array}$$

(19)

where $S_u$ is a function of ${\sigma }_{vi }^{\prime }$ and $OCR$ with means and standard deviations, ${\mu }_{{\sigma }_{vi }^{\prime }}$ and ${\sigma }_{{\sigma }_{vi }^{\prime }}$, ${\mu }_{OCR}$ and ${\sigma }_{OCR}$, respectively, and where $\epsilon$ is the model factor that takes into account the uncertainty of the intrinsic structure model and the effect of sample differences.

2.5. Quantitative Evaluation Indicators

$R^2$, $\mathrm{Evar}$, $\mathrm{RMSE}$, and $\mathrm{MAE}$, the evaluation indicators used in this study, are commonly used to evaluate the prediction performance of the regression model, where $R^2$ and $\mathrm{Evar}$ explain the variance score of the regression model, and their values take the range of [0,1]; values closer to 1 indicate that the independent variable is more able to explain the variance change in the dependent variable, and smaller values indicate the opposite [47,48,49].

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\stackrel{-}{y}\right)}^2}$$

(20)

$$\mathrm{Evar}=\mathrm{Explain} \mathrm{variance}=1-\frac{\mathrm{Var}\left\{y_i-{\widehat{y}}_i\right\}}{\mathrm{Var}\left\{y_i\right\}}$$

(21)

$\mathrm{RMSE}$ indicates the average value of the square root of the error of the fitted data and the original data corresponding to the sample points; $\mathrm{MAE}$ assesses the degree of closeness between the prediction results and the real dataset. The smaller the value of these two, the better the fitting effect. The statistical indicators are computed as follows:

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{n}}$$

(22)

$$\mathrm{MAE}=\frac{{{\displaystyle \sum }}_{i=1}^n\left|y_i-{\widehat{y}}_i\right|}{n}$$

(23)

where $y_i$ denotes the true value of the undrained shear strength of the clay sample, ${\widehat{y}}_i$ denotes the predicted value of the undrained shear strength of the sample clay, $\stackrel{-}{y}$ denotes the value of the undrained shear strength of the clay sample, and $n$ denotes the number of the clay sample. The coefficient of variation (COV) is commonly used in engineering to perform an analysis of variance, showing the variability of a sample population relative to the sample mean.

$$\mathrm{COV}=\frac{\sigma }{\mu }$$

(24)

where $\sigma$ represents the standard deviation of the undrained shear strength of the sample clay, and $\mu$ represents the mean undrained shear strength of the sample clay.

3. Results and Discussion

A comprehensive study of the use of the CatBoost–Bayesian hybrid model with clay-corrected theoretical equations to estimate clay undrained shear strength values is shown in Figure 3.

3.1. Properties of Clays and the Database

The data used in this study were derived from the FI-CLAY/14/856 dataset in the TC304 database, made up of laboratory test data of clay parameters in 33 regions in Finland [50,51]. In this study, 11 parameters of 202 Finnish clay samples were selected for research, including organic content ($Org$), clay content ($CI$), void ratio ($e$), natural water content ($W$), liquid limit ($LL$), plastic limit ($PL$), effective overburden pressure (${\sigma }_{vi}^{\prime }$), preconsolidation pressure (${\sigma }_{v0}^{\prime }$), overconsolidation ratio ($OCR$), compression index ($C_c$), sensitivity ($S_t$), and undrained shear strength ($S_u$); the statistical results of this data set are as follows (

Table 1. Statistical analysis of dataset.

Parameters	Symbol	Min	Max	Std	Mean	COV	Unit
Organic content	$Org$	0.00	7.10	1.48	1.23	1.20	%
Clay content	$CI$	12.70	95.00	20.61	58.99	0.35	%
Void ratio	$e$	0.81	3.88	0.69	2.13	0.32	-
Natural water content	$W$	28.00	155.00	25.52	77.45	0.33	%
Liquid limit	$LL$	24.40	166.00	24.73	76.78	0.32	%
Plastic limit	$PL$	17.70	42.00	4.90	27.68	0.18	%
Effective in situ stress	${\sigma }_{vi}^{\prime }$	4.00	130.00	27.31	41.68	0.66	kPa
Preconsolidation pressure	${\sigma }_{v0}^{\prime }$	13.00	198.00	37.89	69.35	0.55	kPa
Overconsolidation ratio	$OCR$	0.46	20.00	1.81	2.12	0.85	-
Compression index	$C_c$	0.10	4.22	0.86	1.29	0.67	-
Sensitivity	$S_t$	1.69	163	19.97	24.26	0.84	-
Undrained shear strength	$S_u$	5.21	240.00	31.32	28.95	1.08	kPa

Std: Standard deviation.

[51], Figure 4).

Figure 5 shows the normal distribution and fitting effect of organic content ($Org$), clay content ($CI$), void ratio ($e$), natural water content ($W$), liquid limit ($LL$), plastic limit ($PL$), effective overburden pressure (${\sigma }_{vi}^{\prime }$), preconsolidation pressure (${\sigma }_{v0}^{\prime }$), overconsolidation ratio ($OCR$), compression index ($C_c$), and undrained shear strength ($S_u$), and an analysis of the abnormal value of clay undrained shear strength.

Considering the abnormal value of clay undrained shear strength to reduce its influence on the fitting, a simple linear fitting can be performed to obtain the empirical equation of clay undrained shear strength parameters (Figure 5). It could be observed that between the single parameters of clay, the linear relationship is obvious, but the fitting error is large. At the same time, the empirical equation could not be verified. There is a certain particularity, and the constitutive relationship of the clay could not be obtained. Therefore, further research on the relationship between the multi-parameters of the clay is needed.

3.2. Feature Importance of CatBoost–Bayesian Hybrid Model

The CatBoost algorithm (CatBoost–Bayesian hybrid model) was optimized using the Bayesian optimization algorithm under the five-fold cross-validation, and the optimal hyperparameters of the CatBoost algorithm were obtained (

Table 2. CatBoost optimal hyperparameters.

Optuna_Parameters of CatBoost	Description
loss function	MAE
n_estimators	1000
learning rate	0.153
random state	2019
l2_leaf_reg	0.030
colsample_bylevel	0.098
depth	1
boosting type	Plain
bootstrap type	MVS
min_data_in_leaf	4
one_hot_max_size	3
early_stopping_rounds	100

); this study was carried out under the optimal parameters to avoid the excessive model fit, and was beneficial to the CatBoost–Bayesian hybrid model’s adaptively coupled modified theoretical equations for stable estimation of the undrained shear strength of clay. Through the CatBoost–Bayesian hybrid model, an importance ranking of soil parameters that affect the undrained shear strength of clay was obtained.

When the CatBoost–Bayesian hybrid model was used to estimate the undrained shear strength of clay, the feature importance of the model input parameters to the hybrid model under five-fold cross-validation was obtained, and it was explained whether the input parameters contribute positively (positive correlation) or negatively (negative correlation). Each point in the graph represents a data point from the training set. The color represents the value of the feature parameter; red represents the sample with a higher value of the feature parameter, and blue represents the sample with a lower value of the feature parameter. The length of the horizontal line represents the importance of the estimation of the undrained shear strength of the clay.

It can be seen from the whole that the overconsolidation ratio ($OCR$) is the most important parameter for estimating the undrained shear strength of clay, followed by preconsolidation pressure (${\sigma }_{v0}^{\prime }$), effective overburden pressure (${\sigma }_{vi}^{\prime }$), etc. The effective overburden pressure (${\sigma }_{vi}^{\prime }$) increases and the undrained shear strength depends on the consolidation stress before shearing, that is, the research status of the preconsolidation pressure (${\sigma }_{v0}^{\prime }$) [52,53,54]. The void ratio ($e$) and clay content ($CI$) features are the least important. At the same time, the characteristic variables significantly show that the compression index ($C_c$) makes a negative contribution to the computed value of the drainage shear strength, while the preconsolidation pressure (${\sigma }_{v0}^{\prime }$) and the overconsolidation ratio ($OCR$) make a positive contribution to the computed value of the drainage shear strength (Figure 6). The model’s interpretability method is expected to help geotechnical engineers in the selection of soil parameters in practical engineering. However, the CatBoost–Bayesian hybrid model is still a black-box model, which finds it difficult to explain the internal mechanism of the model and the feature importance; therefore, the combination of the feature importance results with the theoretical equations derived from the constitutive model can provide a reference for geotechnical analysis.

3.3. Estimation of Clay Undrained Shear Strength

3.3.1. Uncertainty Analysis of Equation Parameters

Under the Bayesian framework, the CatBoost algorithm (CatBoost–Bayesian) based on a Bayesian optimization algorithm was developed to obtain the feature importance level of soil parameters affecting the undrained shear strength of clay, so as to adaptively couple the theoretical equation of undrained shear strength of $K_0$ consolidated clay, which was derived from the modified Cambridge model; then, the theoretical equation of undrained shear strength of isotropically consolidated clay was established from the critical state of the overconsolidation ratio ($OCR$) and overburden effective pressure (${\sigma }_{vi}^{\prime }$).

Considering the constitutive relation of clay, the calculation results of the established undrained shear strength theoretical equation of isotropically overconsolidated clay were used for inversion calculation. The equation parameter $\gamma$ of the clay in the overconsolidated state is affected by $\left(1-C_s/C_c\right)$, $\gamma =m=0.88\cdot \left(1-C_s/C_c\right)\pm 0.06\cdot SD$, and the equation constant parameter is affected by the plasticity index, $\alpha =2+0.5\cdot I_p$. For the measured parameters of clay that are difficult to obtain, it is recommended that when the clay depth is $1.9 \mathrm{m}<D<17 \mathrm{m}$, the calculation parameters of undrained shear strength should be computed according to the following recommended values ${\left(\alpha \right)}_{\mathrm{OCI}}=2.246$, ${\left(\beta \right)}_{\mathrm{OCI}}=0.490$; the value range of the correction coefficient of ${\sigma }_{vi}^{\prime }$ is $0.441<{\left(\beta \right)}_{\mathrm{OCI}}<0.539$, and the average value is ${\left(\beta \right)}_{\mathrm{OCI}}=0.490$. The theoretical equation is ${\left(S_u\right)}_{\mathrm{OCI}}=2.246\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.490}\cdot {\left(OCR\right)}^{0.770}$. This results in a stable estimate of the undrained shear strength of isotropically consolidated clays (Figure 7a,d).

$${\left(S_u\right)}_{\mathrm{OCI}}=2.246\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.490}\cdot {\left(OCR\right)}^{0.770}$$

(25)

In the normally consolidated state, the fitting parameter ranges were $0.988<{\left(\alpha \right)}_{\mathrm{NCI}}<1.910$, $0.619<{\left(\beta \right)}_{\mathrm{NCI}}<0.786$, and $-0.993<{\left(\gamma \right)}_{\mathrm{NCI}}<0.309$ in this study (Figure 7d). There is a relatively stable estimation result when the clay depth is $1.9m<D<22m$. Therefore, for the calculation of the parameters of undrained shear strength, it is recommended to take the means ${\left(\alpha \right)}_{\mathrm{NCI}}=1.376$, ${\left(\beta \right)}_{\mathrm{NCI}}=0.703$ and ${\left(\gamma \right)}_{\mathrm{NCI}}=-0.342$ (Figure 6b).

$${\left(S_u\right)}_{\mathrm{NCI}}=1.376\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.703}\cdot {\left(OCR\right)}^{-0.342}$$

(26)

In the underconsolidated state, the parameter ranges $0.895<{\left(\alpha \right)}_{\mathrm{UCI}}<3.589$, $0.386<{\left(\beta \right)}_{\mathrm{UCI}}<0.777$, and $0.386<{\left(\beta \right)}_{\mathrm{UCI}}<0.777$ were fitted in this study (Figure 7d). When the clay depth was $2.5 \mathrm{m}<D<17 \mathrm{m}$, the undrained shear strength has a relatively stable estimation result. It is recommended to take the means ${\left(\alpha \right)}_{\mathrm{UCI}}=1.793$, ${\left(\beta \right)}_{\mathrm{NCI}}=0.581$ and ${\left(\gamma \right)}_{\mathrm{NCI}}=-0.423$ for calculation (Figure 7c). The underconsolidated state of clay is the result of the interaction between the strength characteristics of the soil itself and excess pore water pressure, and its undrained shear resistance exhibits structural properties [55]. Therefore, it cannot be verified by the previous results, and it needs to be verified by the relationship between the measured value and the estimated value.

$${\left(S_u\right)}_{\mathrm{UCI}}=1.793\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.581}\cdot {\left(OCR\right)}^{-0.423}$$

(27)

where UCI is the underconsolidated state of isotropic clay.

It can be seen from the whole that the estimated undrained shear strength of clay is affected by the clay depth ($D$) in different consolidation states, which is in line with the actual situation to a certain extent. Asaoka, Guo, and Jiang et al. verified that the undrained shear strength of clay exhibits an obvious linear trend along the depth, through cross-plate shear tests [56,57].

3.3.2. Verification of the Feasibility of the Theoretical Equation

The Bayesian framework concept is an effective means of correcting the original judgments using the new information collected. Based on the new information of the $S_u$ and $OCR$ ofthe Finnish clay database, as well as the ${\sigma }_{vi}^{\prime }$ test index, priori equations are corrected so that the generated posterior equations are more realistic and have fewer errors. The mean and standard deviation of $S_u$ are also calculated to verify the distribution form of the posterior equation.

The form of the probability distribution of $S_u$ for the posterior equation was performed to verify the rationality of the Bayesian framework. By testing, the posterior equation estimated the undrained shear strength with small standard deviation and low variability (Figure 8).

The results visualized in Figure 9 were computed by the theoretical equation mentioned above. The measured and theoretically computed values of the undrained shear strength of clay were compared, in which the abscissa is the measured value, and the ordinate is the theoretical computed value. For the underconsolidated state of clay (Figure 9a), since there was no previous result to verify, Equation (10) of Ohta and Wang et al. was selected to better estimate the undrained shear strength of clay. The theoretically computed value is in good agreement with the measured value. For the normally consolidated state of clay (Figure 9b), the theoretically computed results in this study are slightly larger or smaller than the measured values; meanwhile, for the overconsolidated state of clay (Figure 9c), the theoretically computed results in this study are in good agreement with the measured values, and the computed values of Ohta and Wang et al. are obviously larger. It is found that the $R^2$ values in the underconsolidated state, the normally consolidated state and the overconsolidated state are 0.88, 0.91 and 0.97, respectively.

Under the double logarithmic coordinate, the undrained strength ratio $S_u/{\left({\sigma }_{vi}^{\prime }\right)}^{\beta }$, corrected by $\beta ,$ and the overconsolidation ratio $OCR$ both show an obvious linear relationship, which verifies the theory proposed by our predecessors. $S_u/{\left({\sigma }_{vi}^{\prime }\right)}^{\beta }$ increases linearly with the increase in $OCR$ in the overconsolidated state (Figure 10a), and $S_u/{\left({\sigma }_{vi}^{\prime }\right)}^{\beta }$ decreases linearly with the increase in $OCR$ in the underconsolidated state and normally consolidated state (Figure 10b). The reliability of the theoretical equation proposed in this study is further verified. It is not difficult to find that the gap between the theoretical computed value and the experimental point in this study is significantly smaller than the gap between the theoretically computed value and the measured value studied by Ohta and Wang et al. (Figure 10a).

This study found that ${\left(S_u\right)}_{\mathrm{NCI}}$ has a stable estimation result when the clay depth is $1.9 \mathrm{m}<D<22 \mathrm{m}$, and ${\left(S_u\right)}_{\mathrm{OCI}}$ has stable estimation result when the clay depth is $1.9 \mathrm{m}<D<17 \mathrm{m}$. When the clay depth is $0 \mathrm{m}<D<1.9 \mathrm{m}$, the estimated value of ${\left(S_u\right)}_{\mathrm{OCI}}$ fluctuates greatly. The smaller the compressibility index ($C_c$) value, the lower the compressibility of the soil, and the $C_c$ value of the low compressibility soil is generally less than 0.2. When the depth is less than 1.9 m, the clay in this study is in a state of low compressibility; at the same time, the natural water content is low and the degree of looseness is high, so it is not conducive to estimating the undrained shear strength of the clay. It may be seen from the feature importance ranking of the CatBoost–Bayesian hybrid model that the natural water content ($W$) and the compression index ($C_c$) are important parameters of clay, which can explain the large error in the estimation of ${\left(S_u\right)}_{\mathrm{OCI}}$ by the theoretical equation (Figure 11). Combined with the analysis of the geological tectonic environment, it is because the deposition time of the upper part of the clay depth is relatively short, which is affected by long-term evaporation and water loss. The lower clay has a long deposition time and is affected by the rise and fall in the groundwater level for a long time, which is equivalent to continuous loading and unloading, and finally shows the overconsolidation characteristics and the abrupt change of the clay properties. Therefore, the computed value of the theoretical equation should be much smaller than the measured value.

3.3.3. Comparative Analysis of Estimation Results

Considering the effect of clay depth ($D$), the combination of the feature importance of the CatBoost–Bayesian hybrid model and the calculation of the theoretical equation derived from the modified Cambridge model can well estimate the undrained shear strength of clay, and the theoretical computed value is in good agreement with the measured value (Figure 12a,b). The computed average properties of the clay at different consolidation states were ${\mathrm{R}}^2=0.92$, $\mathrm{Evar}=0.92$, $\mathrm{RMSE}=0.19$, $\mathrm{MAE}=0.03$. This indicated that the uncertainty of the theoretical equation is significantly reduced when the overconsolidation ratio was combined with the preconsolidation pressure or the effective overburden pressure.

The computed results of the test set of the CatBoost–Bayesian hybrid model under five-fold cross-validation are within the 95% confidence interval, and the computed average performance is training set $R^2=0.91$, $\mathrm{Evar}=0.91$, $\mathrm{RMSE}=0.30$, $\mathrm{MAE}=0.12$; test set $R^2=0.86$, $\mathrm{Evar}=0.86$, $\mathrm{RMSE}=0.37$, $\mathrm{MAE}=0.20$. Only a very small number of undrained shear strength estimates at K3 in the test set exceed the 95% confidence interval, proving the reliability of the computed results [58]. In the five-fold cross-validation, there is a gap between the computed value of a small amount of undrained shear strength and the real value, which leads to the large value of $\mathrm{RMSE}$ and $\mathrm{MAE}$, and the fitting effect is not significant enough (Figure 13a–d).

As shown in Figure 14, UCI_$S_u$, NCI_$S_u$, and OCI_$S_u$ demonstrate the performance of CatBoost–Bayesian hybrid model for estimating the undrained shear strength of clay in different consolidation states of clay, as well as the overall CatBoost–Bayesian hybrid model performance, which is subsequently compared with LightGBM-Bayesian, XGBoost-Bayesian hybrid model for comparison. The training and test set performance evaluations are shown in

Table 3. Estimation performance comparison of four hybrid models.

	R2_Mean		Evar_Mean		RMSE_Mean		MAE_Mean
	Train	Test	Train	Test	Train	Test	Train	Test
UCI_Su	-	0.88	-	0.88	-	0.19	-	0.04
NCI_Su	-	0.90	-	0.90	-	0.16	-	0.03
OCI_Su	-	0.97	-	0.97	-	0.22	-	0.02
CatBoost–Bayesian	0.91	0.86	0.91	0.86	0.30	0.37	0.12	0.20
LightGBM-Bayesian	0.99	0.81	0.99	0.81	0.14	0.48	0.08	0.23
XGBoost-Bayesian	0.94	0.81	0.94	0.82	0.24	0.41	0.16	0.24
RandomForest -Bayesian	0.95	0.80	0.95	0.80	0.99	0.98	0.46	0.46

.

${\mathrm{On} \mathrm{the} \mathrm{R}}^2$ and $\mathrm{Evar}$ curves explaining the variance score of the model, the CatBoost–Bayesian hybrid model has the strongest ability to explain the variance in the undrained shear strength under five-fold cross-validation, which could explain almost 86% of the 202 undrained shear strengths of clay. The sample variability shows that the hybrid model has a better effect; $\mathrm{RMSE}$ and $\mathrm{MAE}$ evaluate the closeness of the predicted clay undrained shear strength to its experimental value, and the smaller the value, the better the model fitting effect. It could be observed that the estimated results for clays in different isotropically consolidated states could explain almost 92% of the 202 sample variability in the undrained shear strength of clays. The $\mathrm{RMSE}$ and $\mathrm{MAE}$ of this study are all lower than other intelligent mixed models, and the results show that the combination of CatBoost–Bayesian feature importance and theoretical formula has the best fitting effect (Figure 14).

The estimated performance results in Table 3 show that the theoretical equations of UCI_Su, NCI_Su and OCI_Su derived by combining the CatBoost–Bayesian importance parameters are better than other models in different consolidation states of clay. Among the models of the same type, the estimation performance of the CatBoost–Bayesian hybrid model in the training set and test set is better than other models.

4. Conclusions

Unlike the extensive traditional study of transformation models based on empirical evidence, data-driven ensemble learning methods combined with traditional empirical models have received limited research attention in geotechnical engineering. Moreover, the comprehensive use of Bayesian theory for probabilistic characterization of soil parameter uncertainties and algorithmic optimization problems has not been fully leveraged in a systematic and coherent manner. To address the above challenges, the following study has been conducted. Under the Bayesian framework, the CatBoost algorithm (CatBoost–Bayesian) based on Bayesian optimization algorithm was developed to obtain the feature importance level of soil parameters affecting the undrained shear strength of clay, so as to adaptively couple the theoretical equation of undrained shear strength of $K_0$ consolidated clay, which was derived from the modified Cambridge model, and then the theoretical equation of undrained shear strength of isotropically consolidated clay was established from the critical state of clay parameters, and the calculation results were verified later.

1.

From the feature importance ranking of the CatBoost–Bayesian hybrid model, parameters with high importance and ease of measurement were selected; the overconsolidation ratio ($OCR$) and the effective overburden pressure (${\sigma }_{vi}^{\prime }$) could reasonably explain the model and indirectly estimate the undrained shear strength of the clay.

2.

The equation parameter $\gamma$ of the clay in the overconsolidated state was affected by $\left(1-C_s/C_c\right),\gamma =m=0.88\cdot \left(1-C_s/C_c\right)\pm 0.06\cdot SD$, and the equation parameter $\alpha$ was affected by the plasticity index, $\alpha =2+0.5\cdot I_p$. For the measured parameters of clay that were difficult to obtain, it was recommended that when the clay depth is $1.9m<D<17m$, the calculation parameters of undrained shear strength should be computed according to the following recommended values ${\left(\alpha \right)}_{\mathrm{OCI}}=2.246$, ${\left(\beta \right)}_{\mathrm{OCI}}=0.490$ and ${\left(\gamma \right)}_{\mathrm{OCI}}=0.770$. The theoretical equation was ${\left(S_u\right)}_{\mathrm{OCI}}=2.246\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.490}\cdot {\left(OCR\right)}^{0.770}$.

3.

When the undrained shear strength of clay in the normally consolidated state was estimated at a depth of $1.9m<D<22m$, the recommended theoretical equation was ${\left(S_u\right)}_{\mathrm{NCI}}=1.376\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.703}\cdot {\left(OCR\right)}^{-0.342}$; when the undrained shear strength of clay in the underconsolidated state was estimated at a depth of $1.9m<D<22m$, the recommended theoretical equation was ${\left(S_u\right)}_{\mathrm{UCI}}=1.793\cdot {\left({\sigma }_{vi}^{\prime }\right)}^{0.581}\cdot {\left(OCR\right)}^{-0.423}$.

4.

Compared with the calculation results of Ohta and Wang et al., it was found that the theoretical equation in this study can well estimate the undrained shear strength of isotropically consolidated clay. When the clay depth is $0 m<D<1.9 m$, the huge fluctuation of the estimated value of ${\left(S_u\right)}_{\mathrm{OCI}}$ is mainly due to the long-term influence of evaporative water loss in the upper part of the clay.

5.

The CatBoost–Bayesian hybrid model could excavate the intrinsic relationship of the soil parameters, but it could not give a comprehensive interpretability. The undrained shear strength of isotropic clays was estimated and is to a certain extent interpretable by the CatBoost–Bayesian hybrid model feature importance, adaptively coupled to the theoretical equation derived from the modified Cambridge model. Comparing the results of the CatBoost–Bayesian hybrid model and its similar hybrid models, this study ensured that the average $R^2$ reaches 0.92, the average $\mathrm{RMSE}$ and $\mathrm{MAE}$ were 0.19 and 0.03, respectively, and the overall performance was good.