Prediction of skirted foundation safety factors under combined loading in spatially variable soils using machine learning

Abstract

Skirted foundations are usually used in marine engineering. More researches revealed that the variations in soil undrained shear strength considerably influence the assessing performance of the bearing capacity of skirted foundations. This study proposes two machine learning-based methods to predict safety factors (F_{\mathrm{s}}) of skirted foundations under combined loadings. By comparing the prediction performance of models based on Convolutional Neural Networks (CNN) and Gaussian Process Regression, this study investigates the effect of input size of soil random field on prediction accuracy and identifies the optimal CNN model. The proposed CNN model efficiently predicts corresponding safety factors for different combined loadings under various soil random fields, achieving similar accuracy to the traditional time-consuming random finite element. Specifically, the coefficient of correlation exceeds 0.93 and the mean relative error is less than 2.8% for the variation of the horizontal scales of fluctuation under different combined loadings. The relative error of the predicted F_{\mathrm{s}} value is less than 3.00% given three failure probabilities considering the variation of the vertical scales of fluctuations. These results demonstrate satisfactory prediction performance of the proposed CNN model.

Graphical abstract

Keywords

spatial variability / skirted foundation / safety factor / convolutional neural network / gaussian process regression

1 Introduction

Offshore energy, recognized as a clean and renewable source, has increasingly become a major focus in numerous regions. Because of the complex seabed environment, it is normal to represent the bearing capacity of the foundations under combined loadings through failure envelopes [1]. Significant attention has been directed toward failure envelopes of foundations with combined loadings derived from numerical analysis. However, most existing research primarily focuses on various offshore foundations [2–5].

Due to the finite reserves of resources, marine alternatives present a more abundant alternative. As competition for marine resources intensifies, the demand for various offshore engineering projects has significantly increased [6]. Consequently, skirted foundations have been commonly developed in marine engineering, due to the existence of a peripheral skirt enhances the foundation’s capacity to resist tilting compared to typical surface foundations [7,8]. Thus, the probabilistic failure envelope of skirted foundations under different combined loadings has become a focal point in assessing their ultimate bearing capacity. The accurate prediction of the probabilistic failure envelope is critical, and has garnered growing attention from geotechnical engineers [9–12]. Traditionally, the failure envelope is obtained by deterministic methods, which often neglect the spatially variable soils. In fact, soil properties inherently exhibit spatial variability. Therefore, this spatial variability must be considered when examining probabilistic failure envelopes of foundations.

Soil spatial variability is widely modeled using the random field theory [13–16]. Over recent decades, more studies have been conducted in assessing the stability of slopes [17–20], the stability of tunnel [21–23], and particularly the safety of offshore foundation. For example, Additionally, Fan et al. [24] employed a convex polynomial optimization strategy to represent failure envelopes in spatially variable soils. However, probabilistic failure envelopes derived from traditional methods often require substantial computational effort and are typically limited to obtaining the safety factor (F_{\mathrm{s}}) under specific conditions of spatially variability in soils. These methods generally lack the capacity to effectively analyze the safety performance of skirted foundations under various scenarios.

Owing to the constraints of traditional methods, machine learning has gained more attention and been widely used to solve geotechnical problems, thanks to its development and ability to model nonlinear and extremely complex functions [25–32]. As the probabilistic failure envelope of foundations has increasingly attracted the attention of researchers, the integration of machine learning with traditional research methods, such as field tests, theoretical derivation, and finite element method (FEM), has emerged as a research trend. For example, Chen et al. [33] integrated finite element limit analysis (FELA) with a gradient boosting machine to provide a more accurate and conservative method for predicting failure envelopes. Beygi et al. [34] predicted the bearing capacity of footing of skirted strip by FELA-deep neural network framework. Currently, research primarily focuses on assessing the uniaxial bearing capacity of skirted foundations, with the prediction of F_{\mathrm{s}} corresponding to the failure envelopes of foundations remaining largely unexplored.

This study proposes the application of Convolutional Neural Networks (CNN) and Gaussian process regression (GPR) models to predict the F_{\mathrm{s}} of skirted foundations under combined loadings in spatially variable soils. Spatial variability of soil parameters is generated using random field theory for numerical simulations of the skirted foundation. Then, the failure envelope and corresponding F_{\mathrm{s}} are obtained to create a data set for the subsequent models. Two deep learning models based on CNN and GPR are developed to predict the F_{\mathrm{s}} under different soil random fields. Finally, the effects of the input size and the scales of fluctuation of the random field on model’s prediction performance are estimated.

2 Methodology

2.1 Random field theory

Soil properties typically exhibit spatial variability due to the complex decomposition history [35,36]. The traditional simplified strategy, which assumes the soil as homogeneous material, often leads to discrepancies between calculated results and actual outcomes in the numerical simulation of practical engineering projects. Since Vanmarcke [37] first introduced the random field theory, the simulation of soil spatial variability has increasingly been adopted in the literature for solving geotechnical problems.

It is essential to predefine the relevant location to evaluate the spatial correlation between points within a random field. Since the undrained shear strength (s_{\mathrm{u}}), is widely recognized as significant factor affecting the ultimate bearing capacity of foundations, assumed as the spatially variable parameter. This parameter is modeled using random field theory, as described by Eq. (1):

(1)\rho =\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{2\left|{\tau }_{\mathrm{h}}\right|}{L_{\mathrm{h}}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{2\left|{\tau }_{\mathrm{v}}\right|}{L_{\mathrm{v}}}\right),

where \rho is the correlation function, {\tau }_{\mathrm{h}} and {\tau }_{\mathrm{v}} are the lag distance in the horizontal and vertical directions, L_{\mathrm{h}} and L_{\mathrm{v}} are the autocorrelation distances in the horizontal and vertical directions.

The simulated s_{\mathrm{u}} of soft clay is modeled as adhering to lognormal distribution, aligning with findings that s_{\mathrm{u}} for most submarine clays exhibits either normal or lognormal distributions in reality [38]. The spectral representation method (SRM) is used for generating random fields. The mean of s_{\mathrm{u}} is varied in depth, as expressed by Eq. (2):

(2){\mu }_{s_{\mathrm{u}}}=k\mathrm{z}+{\mu }_{s_{\mathrm{u}0}},

where z is the depth of the soil, {\mu }_{s_{\mathrm{u}}}is the mean value of s_u at the depth z, {\mu }_{s_{\mathrm{u}0}} is the mean value of s_{\mathrm{u}} and set to be 10 kPa, k is the increasing gradient of strength and is indicate of the degree of non-stationarity of this random field, fixed at 1 kPa/m.

The standard Gaussian random field is transformed into the desired non-Gaussian random field using the mean and coefficient of variation (COV), as described in Eqs. (2) and (3):

(3)y\left(x,\mathrm{z}\right)=\overline{y}\left(x,\mathrm{z}\right){\mu }_{s_{\mathrm{u}}}\cdot \mathrm{C}\mathrm{O}\mathrm{V}+{\mu }_{s_{\mathrm{u}}},

where y\left(x,\mathrm{z}\right) is the target non-Gaussian random field, \overline{y}\left(x,\mathrm{z}\right) is the normalized non-Gaussian random field, the COV is set at 0.3.

Phoon and Kulhawy [39] observed that the horizontal scale of fluctuation (SOF_h) is in the range of 40–60 m (with a mean value from the reported literature of 50.7 m). The vertical scale of fluctuation (SOF_v) is in the range of 2–6 m (with a mean value of 3.8 m). Considering the effect of SOFs on the probabilistic failure envelopes of foundations, nine types of random fields with varying SOF_h and SOF_v are simulated. The details of these random filed are listed in Table 1.

Figure 1 demonstrates a typical series of the random field of s_{\mathrm{u}} under various combinations of SOF_h and SOF_v. In the case of unidirectional variation of the random field, it becomes less fluctuating and more stable as the value increases.

2.2 Numerical model

Butterfield et al. [40] proposed a sign convention in which the vertical, horizontal and moment loads were denoted as V, H, and M, respectively. Correspondingly, the vectors of these loadings were presented as u, ω, and θ (as shown in Fig. 2). The subscript “det” or “ran” is used to distinguish results obtained from deterministic or random cases. In this study, the bearing capacity factors are described in Eqs. (4)–(6):

(4)N_{\mathrm{c}\mathrm{H}}=\frac{H_0}{D{\mu }_{{\mathrm{s}}_{\mathrm{u}0}}},

(5)N_{\mathrm{c}\mathrm{V}}=\frac{V_0}{D{\mu }_{{\mathrm{s}}_{\mathrm{u}0}}},

(6)N_{\mathrm{c}\mathrm{M}}=\frac{M_0}{D^2{\mu }_{{\mathrm{s}}_{\mathrm{u}0}}},

where N_cH, N_cV, and N_cM are the bearing capacity coefficient of horizontal, vertical, and moment loading. H₀, V₀, and M₀ are the bearing capacity of pure horizontal, vertical, and moment loading.

The FEM was utilized to analyze the failure envelopes of foundations under small displacement plane strain conditions. The FEM was constructed using version 6.14 of ABAQUS, containing a skirted foundation and surrounding soils. The skirted foundation was modeled as discrete rigid body, with a diameter of 10 m and the thickness of 0.003 m. A reference point (RP) was established at the center of the foundation to apply forces, as shown in Fig. 2. Additionally, the foundation was installed in the soil, with the soil-foundation interface constrained using a “tie” condition [41].

The soil domain was represented with dimensions of 70 m × 50 m to minimize the influence of the boundaries. Horizontal constraints were applied to the lateral boundaries, while both vertical and horizontal constraints were imposed on the bottom boundary. The top boundary remained unconstrained. To enhance model’s computational accuracy, different regions were meshed accordingly. Four-node bilinear reduced integration quadrilateral elements were used for discretization, with a grid size of 1.0 m × 1.0 m applied to the soil region. However, in an area extending 3.0 m × 2.0 m near the foundation, the mesh size was refined to 0.5 m × 0.5 m to mitigate high-stress concentrations at the foundation edges to improve the accuracy of the FEM. In total, the model consisted of 5300 elements, as shown in Fig. 3.

The ideal plasticity model of linear elasticity was applied to the soil, which was treated as an elastic-plastic Mohr–Coulomb material. The Poisson’s ratio of the soil was assumed to be 0.49 and Young’s Modulus was assumed to be E=500s_{\mathrm{u}} [42].

In the deterministic case, the combined loads (u/D=0.07, \omega /D=0.07, \theta /D=0.007) were applied to the RP of the skirted foundation. The calculated N_{\mathrm{c}\mathrm{H},\mathrm{d}\mathrm{e}\mathrm{t}} and N_{\mathrm{c}\mathrm{M},\mathrm{d}\mathrm{e}\mathrm{t}} are 4.20 and 2.54, respectively. The relative errors are less than 3% compared with the results of Selmi et al. [43], (N_{\mathrm{c}\mathrm{H},\mathrm{d}\mathrm{e}\mathrm{t}}=4.07, N_{\mathrm{c}\mathrm{M},\mathrm{d}\mathrm{e}\mathrm{t}}=2.49), indicating that the proposed numerical model demonstrates acceptable accuracy despite slight differences in model mesh sizes.

The displacement-controlled probing method proposed by Bransby and Randolph [44] applies a constant load ratio, beginning at zero load and continuing until load no longer changes with increasing displacement to define the failure envelopes under combined loadings. Constant ratio values of \theta D/\omega, \theta D/u, and u/\omega (unidirectional load is selected as the ultimate load value) are applied at the RP until the ultimate load capacity of the foundation does not change with the increase of displacement, thus obtaining the failure envelopes of the foundations under HM combined loading (where vertical force was not applied), VM combined loading (where horizontal force was not applied), and VH combined loading (where moment was not applied).

Figure 4 presents the failure envelope of the skirted foundation obtained by FEM under HM combined loading excluding the random case analysis. The envelope’s shape closely aligns with the findings of Gourvenec and Barnett [3].

In oceanic foundation projects, the F_{\mathrm{s}} is commonly used to evaluate project safety. The probability of failure is defined as the event where the random result is less than or equal to the deterministic ultimate result divided by F_{\mathrm{s}}, as expressed in Eq. (7) [45]

(7)P_{\mathrm{f}}=P\left(N_{\mathrm{c},\mathrm{r}\mathrm{a}\mathrm{n}}\le \frac{N_{\mathrm{c},\mathrm{d}\mathrm{e}\mathrm{t}}}{F_{\mathrm{s}}}\right),

where N_c,det and N_c,ran are bearing capacity coefficient from deterministic or random case. P_f denotes the probability of occurrence of an event that meets a specific requirement.

According to Eq. (7), calculating F_{\mathrm{s}} typically requires a large number of analyses using the random finite element method (RFEM) to account for spatially variable soils. To enhance prediction efficiency, this study proposes a method whereby F_{\mathrm{s}} can be obtained by computing the relation between the failure envelopes in the deterministic case and random case. The maximum value of F_{\mathrm{s}} is then selected, corresponding to all points on the envelope that simultaneously satisfy this condition, as given by Eqs. (8) and (9):

(8)F_{\mathrm{s}\mathrm{m}}=\frac{N_{\mathrm{c}\mathrm{m},\mathrm{d}\mathrm{e}\mathrm{t}}}{N_{\mathrm{c}\mathrm{m},\mathrm{r}\mathrm{a}\mathrm{n}}},

(9)F_{\mathrm{s}\mathrm{n},\mathrm{r}\mathrm{a}\mathrm{n}}=\mathrm{M}\mathrm{A}\mathrm{X}(F_{\mathrm{s}1},F_{\mathrm{s}2},\cdots ,F_{\mathrm{s}\mathrm{m}}),

where F_{\mathrm{s}\mathrm{m}} is the value of the mth point’s value of the failure envelope, in one case. F_sn,ran is the safety factor from random case.

Figure 5 presents the failure envelopes considering different load ratios under HM combined loading. It is noteworthy that the scaled curve is adjusted by F_{\mathrm{s}\mathrm{n},\mathrm{r}\mathrm{a}\mathrm{n}} based on the deterministic curve.

2.3 Random finite element method

The RFEM is employed in this study to address the challenges posed by spatially variable soils and proposes a non-intrusive stochastic FEM approach by integrating Python with the commercial FEM software ABAQUS. A deterministic FEM of the skirted foundation coupled with soil is developed in the numerical model, with each of the 5300 elements uniquely numbered. The spatially variable s_{\mathrm{u}} is generated using the random field theory aforementioned. These random field values are then assigned to the corresponding elements in the FEM according to their unique identification numbers. Subsequently, an appropriate number of Monte Carlo simulations are selected, and batch calculations are performed using Python scripts to efficiently execute the analyses. The stochastic FEM is used to stochastically analyze F_{\mathrm{s}} of the foundation in random cases. The following probabilistic analysis is conducted on F_{\mathrm{s}} under the nine different cases listed in Table 1.

Figure 6 presents the F_{\mathrm{s}} for the skirted foundation under HM combined loading for Ani-No. 1. It is converged for the mean and standard deviation of F_{\mathrm{s}} after approximately 600 simulations, as the variations in both indicators are less than 2% [46]. Eventually, 800 Monte Carlo simulations are adopted in the random analysis to provide a sufficient number of data sets for following deep learning models.

Figure 7 presents the flow chart of this study. The probabilistic failure envelopes under different combined loadings considering spatially variable soils were calculated using RFEM. Subsequently, the F_{\mathrm{s}} of the probabilistic failure envelopes was calculated using Eqs. (8) and (9). Two surrogated models, based on CNN and GPR, were then developed to predict F_{\mathrm{s}}. Finally, the effect of SOFs on the prediction performance of the optimal surrogated model was analyzed.

3 Machine learning-based method

This study introduces machine learning and deep learning-based models, specifically CNN and GPR, to predict the F_{\mathrm{s}} of foundations in spatially variable soils efficiently and accurately. The inputs for these models consist of a random field data matrix mapping to the soil, analogous to the mapping used in RFEM. The output for the models consists of the F_{\mathrm{s}} of skirted foundations, calculated by RFEM.

3.1 Convolutional neural networks model

CNNs are widely used in geotechnical engineering for classification and regression tasks, attributed to their efficiency and performance [47,48]. This study develops a CNN model on the PyTorch platform to predict the F_{\mathrm{s}} of the skirted foundation under combined loadings The architecture includes an input layer, convolutional layers, pooling layers, fully connected layers (FC layers), and an output layer, as shown in Fig. 8. The input for the model consists of a random field data matrix with dimensions 70 × 50, mirroring the format used in RFEM. The data set is split into training, validation, and testing data sets in 80%, 10%, and 10%. For example, in the scenario of Ani-No. 1, the CNN model is trained on 640 samples, with 80 samples used for validation and another 80 for testing.

The CNN model comprises four convolutional layers that are specifically designed to extract feature information from the input data set. In this research, several hyperparameters of the convolutional layers, such as kernel size, stride, padding, and filters, have been optimized through iterative modifications of the CNN models to achieve optimal predictive performance. Once these hyperparameters are optimized, they are fixed for all subsequent simulations. The kernel size, which defines the receptive field of each convolution operation, is established at 2 × 2 for the initial three layers and 1 × 1 for the final layer. The stride, determining the step size between adjacent convolution operations in both horizontal and vertical directions, is consistently set at 1 × 1 across all layers. Zero padding is applied in the convolution layers to ensure consistent input and output dimensions. The filters control the number of output channels during the convolution process. The activation function used is the rectified linear unit (ReLU), defined by Eq. (10):

(10)f(x)=\{\begin{array}{ll}x,& x>0,\\ 0,& x\le 0.\end{array}

Additionally, the pooling layer plays a critical role in reducing the dimensions of the input features. Based on preliminary results, average pooling was selected for this study, with a pooling size and stride of 2 × 2, effectively halving the input size along the directions. The FC layer, typically positioned at the end of the CNN, facilitates the connection between input and output data sets, mapping the extracted features to the target variable. The expression for the output {\mathit{y}}_{\mathrm{f}} of the FC layer is given by Eq. (11):

(11){\mathit{y}}_{\mathrm{f}}=f({\mathit{W}}_{\mathrm{f}}{\mathit{x}}_{\mathrm{f}}+{\mathit{b}}_{\mathrm{f}}),

where {\mathit{x}}_{\mathrm{f}} and {\mathit{y}}_{\mathrm{f}} are the input and output vectors, {\mathit{W}}_{\mathrm{f}} and {\mathit{b}}_{\mathrm{f}} are the weight and bias matrices, f(\cdot ) is the ReLU activation function described.

To prevent overfitting, early stopping is implemented when the CNN model is trained.

3.2 Gaussian process regression model

The GPR model is increasingly applied as a nonparametric supervised machine learning method in engineering due to its robust theoretical underpinnings in statistical learning and its adeptness at managing complex issues, such as high-dimensional data, small sample sizes, and nonlinear relationships [49–52]. The GPR model is formulated as shown in Eq. (12):

(12)y=f(x)+\epsilon ,

where y is the observation, \epsilon is the Gaussian noise, i.e., \epsilon \sim N(0,{\sigma }_y^2), x is the random variables, and f(x) is the probability density function (PDF) with x, further defined by Eqs. (13)–(15):

(13)f(x)\sim GP(m(x),k\left(x,x^{\prime }\right)),

(14)m(x)=E\left[f(x)\right],

(15)k(x,x^{\prime })=E[(f(x)-m(x))(f(x^{\prime })-m(x^{\prime }))],

where m(x) and k\left(x,x^{\prime }\right) are the mean function and covariance function of the Gaussian process, respectively.

This study introduces a GPR model to predict the F_{\mathrm{s}} of the skirted foundations under combined loadings. The input for the models includes a random field data transformed from two-dimensional to one-dimensional mapping corresponding to the soils. For the GPR model, the data set was 90% for training and 10% for prediction. For example, in the training process of the proposed GPR model using Ani-No. 1, the training data set comprises 720 samples, with 80 samples dedicated to the prediction data set.

The Gaussian kernel function, derived from the density function of the Gaussian distribution, influences the similarity between two data points in various ways. The GPR model employs the radial basis function kernel as the kernel function, given as Eq. (16):

(16)k\left(x,x^{\prime }\right)=\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{\Vert x-x^{\prime }\Vert }^2}{2l^2}),

where {\Vert x-x^{\prime }\Vert }^2 represents the Euclidean distance between x and x', l is the length scale parameter that controls the smoothness or scale of variation of the function.

Typically, the optimal parameters of the GPR model are determined by maximizing the log-marginal likelihood. Given that sparse observations are used, the maximum likelihood estimates often exhibit significant variability near the optimal solution [53]. The GPR model is optimized by iteratively modifying the model to achieve the best prediction performance. These parameters are kept constant in subsequent studies.

4 Results and discussion

4.1 The influence of input size

The results generated by the RFEM are used to train the CNN and GPR models, facilitating efficient prediction of the F_{\mathrm{s}} in spatially variable soils. The performance of these models is evaluated using the correlation coefficient (R), mean relative error (MRE), mean squared error (MSE) and root mean squared error (RMSE). Although the numerical model must consider boundary effects and is relatively large, the influence of soil parameters on the F_{\mathrm{s}} of the foundation is considerably smaller compared to the overall dimensions of the model. However, the size of the input data set significantly impacts the efficiency and accuracy of predictions made by machine learning models. Excessive irrelevant data can introduce noise, leading to issues such as poor convergence or overfitting. To enhance prediction efficiency and accuracy, this study reduces the input data set size around the skirted foundation while maintaining the original size of the numerical model, selecting fewer soil parameters for the new input data set.

Figure 9 illustrates the failure patterns of the skirted foundation under varying load ratios in HM combined loading from deterministic case. Due to the complete symmetry of soil parameters, the failure surface exhibits a symmetric distribution. The failure zone is primarily concentrated within a 40 × 20 region centered around the RF of skirted foundation, with the major failure surface even confined to a smaller 30 × 15 region. Moreover, as shown in Fig. 10, the overall shape of the failure surface remains largely unchanged when considering spatially variable soils. This indicates that the formation and propagation of the failure surface are predominantly influenced by the soil characteristics within the central 40 × 20 region, even confined to a smaller 30 × 15 region. Therefore, special attention is given to this area when constructing the input data set for the subsequent machine learning models. The input domain is progressively reduced from the initial 70 × 50 region, consistent with the finite element model. A faster reduction strategy is applied when the input domain is large, whereas a more refined selection approach is adopted when the input domain approximates the extent of the failure region observed under deterministic conditions.

Table 2 presents the prediction performance of the CNN model for the foundation under HM combined loading for Ani-No. 5, detailing the impact of varying input data set sizes on the CNN model’s predictive efficacy. It is noted that the CNN model achieves optimal predictive performance with an input data set size of 486, confirming that the optimization of data set size improves prediction efficiency and accuracy.

Figure 11 presents the prediction performance of the CNN and GPR models with varying input data set sizes under four typical conditions of the foundation in HM combined loading (Ani-No. 1, Ani-No. 6, Ani-No. 7 and Ani-No. 9). The graphs indicate that as the data set size decreases from 3500, the overall trends of R and MRE values remain similar. When the input data set is large, the inclusion of marginal soil parameters with minimal influence on F_{\mathrm{s}} significantly hampers the training of both the CNN and GPR models, resulting in longer convergence times, lower R values, higher MRE values, and poorer predictive performance. As the input data set size decreases, the CNN and GPR models both show improved prediction accuracy, with the CNN model outperforming the GPR model, particularly when the data set size is reduced to 486. However, further reductions in data set size result in diminished predictive performance for both models, as critical soil parameters influencing F_{\mathrm{s}} are not adequately represented in the input data.

Additionally, after the optimal hyperparameters were determined, training a single CNN model required approximately 35 min, while the GPR model took around 15–20 min. Although the CNN model requires a longer training time compared to the GPR model, it still represents a significant improvement in computational efficiency relative to the RFEM simulations. Furthermore, the CNN model demonstrated superior predictive performance over the GPR model, particularly in capturing complex nonlinear patterns. Therefore, subsequent analysis in this study focuses on the CNN model.

4.2 The influence of horizontal scale of fluctuation of the random method

In Subsection 4.1, the values of R, MRE, RMSE and MSE are used as evaluation metrics, indicating that the CNN model with an input data set size of 468 outperforms the GPR model. The RMSE values all below 10% of the F_{\mathrm{s}} generally indicates satisfactory predictive accuracy. Therefore, the values of R and MRE were emphasized, as they are intuitive and widely accepted metrics for assessing trend consistency and relative prediction accuracy.

Figure 12 presents the performance of the CNN models under HM combined loading. The R values range from 0.930 to 0.990, all exceeding 0.9338, and the MRE values range from 1.8% to 2.7%, all below 2.66%. These results indicate the strong predictive accuracy of the developed CNN model for F_{\mathrm{s}}of the foundations.

With the SOF_h increasing, the CNN model exhibits improved predictive performance, attributed to stabilization of fluctuations in the random field. Figure 10 illustrates that as the SOF_h increases from 10 to 50 m, R values continuously rise, reaching a maximum of 0.9825 and a minimum MRE of 1.92% at SOF_h = 50 m. Beyond this point, as SOF_h continues to increase, the R value begins to decline. In the analysis of other combined loadings, Fig. 13(a) shows that the R values range from 0.960 to 0.980, all exceeding 0.9614. The MRE values range from 1.8% to 2.6%, which the MRE values are smaller than 2.53%. Similarly, Fig. 13(b) presents R values are all exceeding 0.9396 and the MRE values are all below 2.42%. There is a consistent pattern where the value of R reaches its maximum at SOF_h = 50 m under three combined loadings.

For the foundation under three combined loadings, further analysis of the predictive performance at SOF_h = 50 m and SOF_h = 60 m is detailed in Fig. 14. The results at SOF_h = 60 m surpass those at SOF_h = 50 m without considering the impact of input data set size. However, when the input data set size is factored in, the results at SOF_h = 50 m excel compared to those at SOF_h = 60 m. The consistent patterns across the three combined loadings suggest that this behavior is systematic and not merely a random occurrence due to the stochastic distribution of soil parameters. Further investigation is warranted into the different distribution patterns of spatially variable soils following modifications in the dimensions of the input data set.

Figure 15 shows the random field of s_{\mathrm{u}} of soil, corresponding to the minimum and maximum relative R values of Ani-No. 5 and Ani-No. 6 under HM combined loading. The outer edge dimensions of Fig. 14 represent the finite element soil region, while the highlighted areas reflect the range of soil parameters after the data set optimization. These random field indicate a more uniform stratified distribution of soil parameters at SOF_h = 50 m compared to SOF_h = 60 m, resulting in more concentrated F_{\mathrm{s}} values. Consequently, the predictive performance is superior for SOF_h = 60 m.

Columns a to d in Table 3 corresponds to the conditions presented in Figs. 15(a)–15(d), and the observed phenomena can be further elucidated from the perspective of the COV. The COV value of soil parameters at SOF_h = 60 m is higher than that at SOF_h = 50 m, indicating greater dispersion in the input soil parameters. This increased variability provides more valuable information for the CNN model during training, thereby enhancing the prediction accuracy for SOF_h = 60 m compared to SOF_h = 50 m. Similarly, when a smaller input data set is considered, the COV value for the SOF_h = 50 m random field is higher, reflecting a greater dispersion of the input soil parameters and thus providing more training information, which leads to better prediction performance.

For scenarios where SOF_h varies from 10 to 60 m, six individual surrogate models with identical structures (as shown in Fig. 8) are developed. These models differ in parameters due to training on separate data sets. However, generating multiple data sets and training separate models is a complex and resource-intensive process, particularly when addressing F_{\mathrm{s}} prediction under varying SOF_h in random fields. A unified surrogate model is thus crucial for stochastic analysis across different SOF_h values, provided it maintains acceptable prediction accuracy. Therefore, the samples from Ani-No. 1 to Ani-No. 6 are combined to create a mixed data set for training a unified CNN model to predict F_{\mathrm{s}} for skirted foundations. The training process involves 3840 samples for model training, 480 samples for validation, and 480 samples for testing.

Table 4 presents the relationship between the CNN models and unified CNN model under different combined loadings. The R values of the unified CNN models exceed the average of CNN models, and exceeds all the values of CNN models considering the variation of SOF_h in the random field under combined loadings. The MRE values of the unified CNN models are smaller than the average of CNN models and all the values of CNN models under different combined loadings. The MRE values of the unified CNN models are smaller than 2.00%.

Figures 16(a), 17(a), and 18(a) illustrate the PDF of the predicted and actual values of F_{\mathrm{s}} under HM combined loading, VM combined loading and VH combined loading. Across all loading combinations, the relative errors remain below 10%. The confidence intervals (CIs) of the relative errors are 95.8%, 95.8% and 95.6%, with a margin of error of 5%. Figures 16(b), 17(b), and 18(b) display the cumulative density function (CDF) of the predicted and actual values of F_{\mathrm{s}} under various combined loadings. The PDFs of the predictions generated by the unified CNN model align closely with the actual values derived from RFEM, demonstrating strong consistency. These confirm that the proposed unified CNN model effectively captures the global distribution of the F_{\mathrm{s}} while accounting for the variability of SOF_h of the random field. Additionally, the mean, minimum, and maximum values of predicted and actual values are compared to evaluate predictive performance. Additionally, comparisons of the minimum, maximum and mean predicted and actual values further evaluate the model’s predictive performance. For HM combined loading, these relative errors are 0.51%, 2.14%, and 0.54%. For VM combined loading, these relative errors are 2.61%, 1.18%, and 0.01%. For VH combined loading, hese relative errors are 2.21%, 1.39%, and 0.44%. All relative errors remain within 3.00%.

Table 5 presents the relationship between the predicted results and actual values, denoted as “Pre” and “Act”, considering the 10%, 50%, and 90% envelopes of unified CNN model under combined loadings. The true results are calculated using RFEM, and the predicted results are predicted through the CNN model. The relative error values range from 0.86% to 1.15%. These results demonstrate that the unified CNN model provides accurate prediction performance in estimating the F_{\mathrm{s}} of foundations in spatially variable soils.

4.3 The influence of vertical scale of fluctuation of the random method

In Subsection 4.2, R and MRE values are employed as evaluation metrics. However, due to the number of testing data sets, MRE does not exhibit a clear variation pattern and only partially reflects the prediction performance of the CNN models. Moreover, it is sensitive to outliers. Therefore, this section primarily uses the R value as the evaluation metric, as it more accurately reflects the predictive performance of the model.

Figure 19 presents the distribution of R values in different SOF_v soils for the foundation under HM combined loading. The R values range from 0.955 to 0.985, all exceeding 0.9598, indicating that the CNN model can predict the F_{\mathrm{s}} of foundation well. As the SOF_v increases, the R values rise, reaching a peak of 0.9843 at SOF_v = 8 m. This trend is attributed to the fact that an increase in SOF_v leads to a decrease in the amplitude of random field fluctuations, resulting in more concentrated soil parameter values and a more uniform distribution of F_{\mathrm{s}}. Consequently, the model is better able to capture the relationship between soil parameters and F_{\mathrm{s}} when predicting.

Figure 20 illustrates the distribution of R values in different SOF_v under VM combined loading. The R values range from 0.965 to 0.990, with the smallest value of 0.9660. As the SOF_v increases, R values also increase, indicating an improvement in the model’s predictive performance.

Similarly, Fig. 21 presents the distribution of R values in different SOF_v under the foundation subjected to the VM combined loading. The R values range again from 0.965 to 0.990, with the smallest value of 0.9668. As the SOF_v increases, there is a corresponding increase in R values, indicating a consistent improvement in prediction accuracy. The model performance exhibits a similar trend to the previous cases, with R values increasing as the SOF_v increases. This consistency underscores the model’s enhanced predictive capability as the SOF_v increases

The F_{\mathrm{s}} is indicative of the bearing capacity of the foundation both in deterministic case and random scenarios. When the F_{\mathrm{s}}>1.0, it indicates that the foundation requires a higher bearing capacity to ensure safety in spatially variable soils.

Table 6 presents the relationship between the predicted and actual values of F_{\mathrm{s}} considering the 10%, 50%, and 90% envelopes in HM combined loading considering the SOF_v. The largest relative error recorded is 2.53%.

Table 7 presents the relationship between the predicted results and actual values of F_{\mathrm{s}} considering the 10%, 50%, and 90% envelopes in VM combined loading considering the SOF_v. The largest relative error recorded is 2.25%.

Table 8 presents the relationship between the predicted and actual values of F_{\mathrm{s}} considering the 10%, 50%, and 90% envelopes in VH combined loading considering the SOF_v. The largest relative error recorded is 2.52%.

The relative errors of the CNN models in predicting F_{\mathrm{s}} are less than 3.00% under three combined loadings. Thus, it can be concluded that the CNN models can effectively replace the time-consuming RFEM model for predicting the F_{\mathrm{s}} of foundations in spatially variable soils.

For scenarios where SOF_v varies from 2 and 8 m, four individual surrogate models with identical structures (as shown in Fig. 8) are developed. These models differ in parameters due to training on separate data sets. However, generating multiple data sets and training separate models is a complex and resource-intensive process, particularly when addressing F_{\mathrm{s}} prediction under varying SOF_v in random fields. A unified surrogate model is thus crucial for stochastic analysis across different SOF_v values, provided it maintains acceptable prediction accuracy. Therefore, the Ani-No. 5, Ani-No. 7, Ani-No. 8 and Ani-No. 9 are combined to create a mixed data set for training a unified CNN model to predict the F_{\mathrm{s}} for skirted foundations. The training process involves 2560 samples for model training, 320 samples for validation, and 320 samples for testing.

Table 9 presents the R and MRE between the CNN models and unified CNN model under different combined loadings. The R values of the unified CNN models exceed the average of CNN models, and exceeds all the values of CNN models considering the variation of SOF_v in the random field under different combined loadings. The MRE values of the unified CNN models are smaller than the average of CNN models and all the values of CNN models under combined loadings. The MRE values of the unified CNN models are smaller than 2.00%.

Table 10 presents the relative errors remain below 10%, and the CIs of the relative errors are 96.3%, 97.2%, and 96.9%, with a margin of error of 5%, across all loading combinations. Additionally, comparisions of the minimum, maximum and mean predicted and actual values further evaluate the model’s predictive performance. For HM combined loading, the relative errors of the maximum, minimum, and mean values between the predicted and actual values are 0.13%, 0.95%, and 0.31%, respectively. For VM combined loading, these relative errors are 0.60%, 1.19%, and 0.68%. For VH combined loading, these relative errors are 1.51%, 1.46%, and 0.10%. All relative errors remain within 2.00%.

Table 11 presents the relationship between the predicted results and actual values of F_{\mathrm{s}} considering the 10%, 50%, and 90% envelopes of unified CNN model under combined loadings. The relative error values range from 0.86% to 1.23%. These results demonstrate that the unified CNN model provides an accurate and efficient alternative to the computationally intensive FEM for estimating the distribution of F_{\mathrm{s}}.

5 Conclusions

In summary, this study integrates random field theory with CNN and GPR to establish an effective machine learning model for predicting the F_{\mathrm{s}} of skirted foundations in spatially variable soils under three combined loadings. The SRM is used to simulate the non-stationary random field of s_{\mathrm{u}} in soil, considering SOF_h and SOF_v, which as input for the CNN and GPR models. The F_{\mathrm{s}} corresponding to the probabilistic failure envelope of combined loadings, is defined as the output of the CNN and GPR models. Comparisons of the predictions with RFEM results validate the model’s performance. Key conclusions include

1) The CNN model can effectively replace the time-consuming RFEM for predicting the F_{\mathrm{s}} of skirted foundations with reasonable accuracy. RFEM requires approximately 1100 h, while the CNN and GPR models only need 4 and 2.5 h to train all models and almost 10 s to predict F_{\mathrm{s}}, respectively. Although the GPR model is faster, the CNN model provides higher accuracy, making it the preferred surrogate model.

2) The CNN model, which uses a small range of soil parameters near the skirted foundation as the input data set containing 468 data points, performs better than models using parameters from the entire area.

3) Three CNN models with the same architecture predict the F_{\mathrm{s}} of skirted foundations under different combined loadings, considering the effect of SOF_h. The CNN model demonstrates the best prediction performance when the input data set size is small, the COV is larger at SOF_h = 50 m, and the input set contains more valid information.

4) Three CNN models predict the F_{\mathrm{s}} considering the effect of SOF_v, which exhibit the best prediction performance at SOF_v = 8 m. The relative errors of the proposed CNN models in predicting F_{\mathrm{s}} are all below 3.00% under three combined loadings in 10%, 50%, and 90% envelope.

5) Two unified CNN models predict the F_{\mathrm{s}} under combined loadings considering the effects of SOF_h and SOF_v. The relative errors of the proposed unified CNN models in predicting F_{\mathrm{s}} are less than 3.00% under three combined loadings.

The proposed CNN model effectively predicts the F_{\mathrm{s}} of skirted foundations under combined loadings in variability spatial soils. However, two important limitations remain and warrant further investigation. First, relatively large prediction deviations may occur at the tails of the F_{\mathrm{s}} distribution, particularly beyond the 95%CI. These deviations are primarily due to the limited representation of low-probability failure events in the training data, which restricts the model’s ability to accurately capture extreme responses. Second, the current model does not address predictions under H–V–M combined loading. Future work will consider data enhancement strategies, such as importance sampling, to improve prediction accuracy in these regions. Besides, the study will focus on developing a generalized deep learning framework capable of handling three-dimensional loading conditions. This will involve defining a new scaling indicator to extend the failure envelope from a curve to a surface, thereby enabling more comprehensive probabilistic predictions.

References

[1] Roscoe K . The stability of short pier foundations in sand. British Welding Journal, 1956, 3(8): 343–354

[2] Bransby M F , Randolph M F . Combined loading of skirted foundations. Geotechnique, 1998, 48(5): 637–655

[3] Gourvenec S , Barnett S . Undrained failure envelope for skirted foundations under general loading. Geotechnique, 2011, 61(3): 263–270

[4] Luo W Z , Li J H . Effects of installation on combined bearing capacity of a spudcan foundation in spatially variable clay. Applied Ocean Research, 2022, 120: 103055

[5] Zou J X , Yang J Z , Hu F J . Bearing capacity of monopile-bucket composite foundation in sand-over-clay under V-H-M combined static loads. Applied Ocean Research, 2024, 150: 104092

[6] Randolph M F , Gaudin C , Gourvenec S M , White D J , Boylan N , Cassidy M J . Recent advances in offshore geotechnics for deep water oil and gas developments. Ocean Engineering, 2011, 38(7): 818–834

[7] Haddad A , Barari A , Amini R . The remedial performance of suction caisson foundations for offshore wind turbines under seismically induced liquefaction in the seabed: Shake table testing. Marine Structures, 2022, 83: 103171

[8] Cheng H J , Zhou M , Niu F J , Zhang X H , Tian Y H , Li J H . Numerical study of the bearing capacity of the hybrid skirted foundations in silty sand under combined loading. Applied Ocean Research, 2024, 150: 104129

[9] Yun G , Bransby M F . The horizontal-moment capacity of embedded foundations in undrained soil. Canadian Geotechnical Journal, 2007, 44(4): 409–424

[10] Bransby M F , Yun G J . The undrained capacity of skirted strip foundations under combined loading. Geotechnique, 2009, 59(2): 115–125

[11] Mehravar M , Harireche O , Faramarzi A . Evaluation of undrained failure envelopes of caisson foundations under combined loading. Applied Ocean Research, 2016, 59: 129–137

[12] Wang Y K , Zou X J , Zhou M , Zhang X H . Capacity and failure mechanism of monopile-wheel hybrid foundation in clay-overlaying-sand deposits under combined V-H-M loadings. Marine Structures, 2023, 90: 103443

[13] Cho S E . Probabilistic assessment of slope stability that considers the spatial variability of soil properties. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(7): 975–984

[14] Huang J , Griffiths D V , Fenton G A . Probabilistic analysis of coupled soil consolidation. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(3): 417–430

[15] Bari M W , Shahin M A . Probabilistic design of ground improvement by vertical drains for soil of spatially variable coefficient of consolidation. Geotextiles and Geomembranes, 2014, 42(1): 1–14

[16] Li K , Wang R , Ma H B , Zhang J M . Rising groundwater table due to restoration projects amplifies earthquake induced liquefaction risk in Beijing. Nature Communications, 2025, 16(1): 1466

[17] Li J H , Tian Y H , Cassidy M J . Failure mechanism and bearing capacity of footings buried at various depths in spatially random soil. Journal of Geotechnical and Geoenvironmental Engineering, 2015, 141(2): 04014099

[18] Jiang S H , Li D Q , Cao Z J , Zhou C B , Phoon K K . Efficient system reliability analysis of slope stability in spatially variable soils using Monte Carlo simulation. Journal of Geotechnical and Geoenvironmental Engineering, 2015, 141(2): 04014096

[19] Li X Y , Liu X , Liu Y D , Yang Z Y , Zhang L M . Probabilistic slope stability analysis considering the non-stationary and spatially variable permeability under rainfall infiltration-redistribution. Bulletin of Engineering Geology and the Environment, 2023, 82(9): 350

[20] Jiang S H , Hu H P , Wang Z Z . Improved Bayesian model updating of geomaterial parameters for slope reliability assessment considering spatial variability. Structural Safety, 2025, 112: 102536

[21] Ali A , Lyamin A V , Huang J S , Sloan S W , Cassidy M J . Undrained stability of a single circular tunnel in spatially variable soil subjected to surcharge loading. Computers and Geotechnics, 2017, 84: 16–27

[22] Chen R P , Zhang P , Wu H , Wang Z , Zhong Z . Prediction of shield tunneling-induced ground settlement using machine learning techniques. Frontiers of Structural and Civil Engineering, 2019, 13(6): 1363–1378

[23] Zhang W G , Han L , Gu X , Wang L , Chen F Y , Liu H L . Tunneling and deep excavations in spatially variable soil and rock masses: A short review. Underground Space, 2022, 7(3): 380–407

[24] Fan S T , Zhang Y R , Li S , Han M L . Probabilistic failure envelopes of tripod bucket-supported offshore wind turbines in spatially variable clay based on a novel continuous method. Computers and Geotechnics, 2024, 170: 106305

[25] Ismail A , Jeng D S , Zhang L L . An optimised product-unit neural network with a novel PSO-BP hybrid training algorithm: Applications to load-deformation analysis of axially loaded piles. Engineering Applications of Artificial Intelligence, 2013, 26(10): 2305–2314

[26] Zhou C , Xu H C , Ding L Y , Wei L C , Zhou Y . Dynamic prediction for attitude and position in shield tunneling: A deep learning method. Automation in Construction, 2019, 105: 102840

[27] Nguyen-Le D H , Tao Q B , Nguyen V H , Abdel-Wahab M , Nguyen-Xuan H . A data-driven approach based on long short-term memory and hidden Markov model for crack propagation prediction. Engineering Fracture Mechanics, 2020, 235: 107085

[28] Deb P , Pal S K . Interaction behavior and load sharing pattern of piled raft using nonlinear regression and LM algorithm-based artificial neural network. Frontiers of Structural and Civil Engineering, 2021, 15(5): 1181–1198

[29] Fu Y P , Lin M S , Zhang Y , Chen G F , Liu Y J . Slope stability analysis based on big data and convolutional neural network. Frontiers of Structural and Civil Engineering, 2022, 16(7): 882–895

[30] Fernández A , Sanchidrian J A , Segarra P , Gomez S , Li E M , Navarro R . Rock mass structural recognition from drill monitoring technology in underground mining using discontinuity index and machine learning techniques. International Journal of Mining Science and Technology, 2023, 33(5): 555–571

[31] Abdollahi A , Li D L , Deng J , Amini A . An explainable artificial-intelligence-aided safety factor prediction of road embankments. Engineering Applications of Artificial Intelligence, 2024, 136: 108854

[32] Kong X R , Wu Y C , Zhu P , Zhi P , Yang Q F . Novel artificial neural network aided structural topology optimization. Applied Sciences, 2024, 14(23): 11416

[33] Chen H Z , Shen Z C , Wang L , Qi C C , Tian Y H . Prediction of undrained failure envelopes of skirted circular foundations using gradient boosting machine algorithm. Ocean Engineering, 2022, 258: 111767

[34] Beygi M , Fallahi M , Vali R , Mousavi E , Saberian M , Li J , Barari A . FELA-DNN framework to predict the seismic bearing capacity of skirted strip footing built on a non-cohesive slope. Soil Dynamics and Earthquake Engineering, 2023, 171: 107932

[35] Wang Y , Cao Z J . Probabilistic characterization of Young’s modulus of soil using equivalent samples. Engineering Geology, 2013, 159: 106–118

[36] Cao Z J , Wang Y , Li D Q . Quantification of prior knowledge in geotechnical site characterization. Engineering Geology, 2016, 203: 107–116

[37] Vanmarcke E H . Probabilistic modeling of soil profiles. Journal of the Geotechnical Engineering Division-ASCE, 1977, 103(11): 1227–1246

[38] Lacasse S , Nadim F . Uncertainties in characteriscing soil properties. In: Uncertainty in the Geologic Environment: From Theory to Practice, 1996, 1(58): 49–75

[39] Phoon K K , Kulhawy F H . Characterization of geotechnical variability. Canadian Geotechnical Journal, 1999, 36(4): 612–624

[40] Butterfield R , Houlsby G T , Gottardi G . Standardized sign conventions and notation for generally loaded foundations. Geotechnique, 1997, 47(5): 1051–1054

[41] Ye Z T , Gao Y F , Shu S , Wu Y X . Probabilistic undrained bearing capacity of skirted foundations under HM combined loading in spatially variable soils. Ocean Engineering, 2021, 219: 108297

[42] Hung L C , Kim S R . Evaluation of combined horizontal-moment bearing capacities of tripod bucket foundations in undrained clay. Ocean Engineering, 2014, 85: 100–109

[43] Selmi M , Kormi T , Hentati A , Bel Hadj A N . Capacity assessment of offshore skirted foundations under HM combined loading using RFEM. Computers and Geotechnics, 2019, 114: 103148

[44] Bransby M F , Randolph M F . The effect of embedment depth on the undrained response of skirted foundations to combined loading. Soil and Foundation, 1999, 39(4): 19–33

[45] Wu Y X , Zhang H L , Shu S . Probabilistic bearing capacity of spudcan foundations under combined loading in spatially variable soils. Ocean Engineering, 2022, 248: 110738

[46] Huang H W , Xiao L , Zhang D M , Zhang J . Influence of spatial variability of soil Young’s modulus on tunnel convergence in soft soils. Engineering Geology, 2017, 228: 357–370

[47] Cha Y J , Choi W , Büyüköztürk O . Deep learning-based crack damage detection using convolutional neural networks. Computer-Aided Civil and Infrastructure Engineering, 2017, 32(5): 361–378

[48] Yin X L , Wang H , Pisanò F , Gavin K , Askarinejad A , Zhou H P . Deep learning-based design model for suction caissons on clay. Ocean Engineering, 2023, 286: 115542

[49] Mahmoodzadeh A , Mohammadi M , Salim S G , Ali H F H , Ibrahim H H , Abdulhamid S N , Nejati H R , Rashidi S . Machine learning techniques to predict rock strength parameters. Rock Mechanics and Rock Engineering, 2022, 55(3): 1721–1741

[50] Tamhidi A , Kuehn N M , Bozorgnia Y . Uncertainty quantification of ground motion time series generated at uninstrumented sites. Earthquake Spectra, 2023, 39(1): 551–576

[51] Ching J , Yoshida I , Phoon K K . Comparison of trend models for geotechnical spatial variability: Sparse bayesian learning vs. gaussian process regression. Gondwana Research, 2023, 123: 174–183

[52] Zhang H L , Wu Y X , Yang S C . Probabilistic analysis of tunnel convergence in spatially variable soil based on Gaussian process regression. Engineering Applications of Artificial Intelligence, 2024, 131: 107840

[53] Li R Z , Sudjianto A . Analysis of computer experiments using penalized likelihood in Gaussian kriging models. Technometrics, 2005, 47(2): 111–120

RIGHTS & PERMISSIONS

Higher Education Press