Prediction and Analysis of Axial Stress of Piles for Piled Raft Due to Adjacent Tunneling Using Explainable AI

Abstract

Tunneling, especially in urban areas, affects many structures on the ground, which directly influences the usability and stability of the structures. The settlement of and axial stress on the pile foundation are important factors that determine the behavioral characteristics of the pile foundation. Therefore, this study uses numerical analysis and machine learning to derive a prediction model of pile axial stress due to tunnel excavation adjacent to the piled raft. Numerical analysis data were utilized for machine learning purposes, and the effects of the input data on the prediction model were scrutinized. The numerical analysis revealed that the change in the pile axial stress resulting from tunnel excavation differed depending on the pile’s location, with the greatest axial stress reduction occurring in the center of the piled raft. Furthermore, the rate of reduction was higher in soils with lower relative densities. Several algorithms were employed to derive the prediction model, with tree-based algorithms displaying notable performance in predicting pile axial stress. Additionally, preprocessing the data with appropriate feature engineering techniques exhibited superior predictive power, and incorporating settlement data aided in enhancing the prediction model’s performance.

1. Introduction

Pile foundations are widely used to support urban structures, and tunnel-induced ground deformation is transmitted to these piles, which in turn transmit the deformation to the superstructure. The variation in pile settlement and axial force has been extensively studied by many researchers as a means of estimating the state of a pile foundation through field monitoring, and the behavioral characteristics of piles affected by tunnel excavation are also being actively investigated [1,2,3,4,5,6,7,8,9].

Since its introduction by Poulos [10], piled raft foundations have been applied to numerous large structures. Although many variables must be considered in piled raft design, traditional regression analysis techniques are currently limited to the consideration of only a few variables.

In a study on load transfer induced by shield tunnelling in soft ground, Zheng et al. [11] investigated changes in the load acting on piles through both field measurements and numerical analysis of the tunnel location. The results of the study indicate that changes in the load due to tunnel excavation are dependent on the pile’s location. Specifically, the load on the pile located in the middle of the foundation decreased when the tunnel was excavated under the pile, while the load on the piles at both ends tended to increase to some extent. The authors of this paper intend to compare the outcomes of this study with the numerical analyses conducted in their research.

Machine learning, with its ability to consider multiple variables simultaneously, is an effective tool for overcoming the limitations of traditional regression analysis techniques and is used in many fields, including engineering, natural sciences, and medicine, as well as geotechnical engineering. Li et al. [12] analyzed single masonry piers with different aspect ratios (h/l) by employing the multi-pier (MP) method. The results of the analysis were used as input to a machine-learning algorithm called back-propagation multilayer perceptron (BPMLP) to generalize the data, and they predicted the potential failure modes of the masonry wall. Rehman et al. [13] presented an intelligent model of the hydraulic conductivity of sandy soil by employing various machine-learning algorithms, among which they proposed a GEP-based model to predict the conductivity value with the best performance. They resolved the limitations of the existing value prediction models by using a broad range of input parameters that define geological characteristics, such as grain size, gradation parameters, and dry density. Onyelowe et al. [14] provided a guide to the application of machine learning to various geotechnical problems by reviewing intelligent learning algorithms such as ANN, Fuzzy Logic, GEP, ANFIS, and ANOVA applied to the prediction of geotechnical problems. Therefore, the performance of the final predictive model is affected by the variables of the input data used, missing values, and hyperparameters of the algorithm, which must be considered when utilizing machine learning. According to Ijaz et al. [15], various geotechnical parameters demonstrate the reliability of developed GSMs through high KPIs (RMSE, MAE, NSE, and PC), and models based on various geotechnical parameters provide efficient design solutions. In particular, according to Andrew Ng [16], the quality of the data is a crucial factor in determining the performance of a predictive model. Khalid et al. [17] proposed several predictive models to estimate the physico-mechanical characteristics of undisturbed soils. Among them, the 3D response surface-analysis-based models showed 8.2–42.4% less deviation for all parameters compared with the model developed using only DCPI as an independent variable.

In the present study, the authors utilized numerical analysis data to develop a model capable of predicting the axial stress induced in piles due to tunnel excavation in the vicinity of a piled raft, through the application of machine-learning techniques. To evaluate the performance of the models, the authors varied the input data in multiple ways and employed various machine-learning algorithms, optimizing the hyper-parameters of each algorithm. In order to overcome the challenge of the “black-box” nature of machine-learning models, which can provide high accuracy at the cost of low interpretability, the authors employed Explainable AI to quantitatively assess the contribution of each input variable to the final prediction model.

2. Numerical Analysis

2.1. Modeling and Material Properties

This study presents a three-dimensional numerical analysis model to simulate a tunnel excavation adjacent to a piled raft and to evaluate the settlement of the raft and the axial stress of the pile. The numerical analysis model encompasses a range of variables, including the modulus of elasticity of the ground, pile diameter, pile length, foundation width, number of piles, pile spacing, pile-to-tunnel horizontal and vertical separation distances, and tunnel face-to-foundation horizontal distance, which are outlined in

Table 1. Variables for numerical analysis.

No.	Variable	Acronym	Unit	Range
1	Relative density	EG	kPa	20 × 10−3, 40 × 10−3, 50 × 10−3
2	Horizontal offset pile–tunnel	OH	m	0.0TD, 2.0TD, 5.0TD (1)
3	Vertical offset pile–tunnel	OV	m	0.0TD, 1.5TD, 5.0TD (1)
4	Pile diameter	PD	m	0.5, 0.8, 1.0
5	Pile length	PL	m	6.5, 10.0, 12.5, 20.0, 22.0, 25.0, 35.0, 44.0
6	Raft width	BR	m	12.0, 20.0, 25.0
7	Pile number	PN	ea	9, 25
8	Pile spacing	PCTC	m	1.25PD, 2.5PD, 5.0PD (2)
9	Offset raft–tunnel face	TF	m	10.0BR, 4.0BR, 1.0BR, 0.0BR (3)

¹ TD: Tunnel diameter; ² PD: Pile diameter; ³ BR: raft width.

along with the scope of their application.

The elastic modulus of the ground is classified based on the relative density of the soil, as presented by Lambe and Whitman [18]. Furthermore, the pile diameter and pile length are categorized as either long piles or short piles, as per the guidelines outlined in the Basic Design Criteria for Structures [19]. The pile spacing is considered both as a group pile and a single pile by referring to the Basic Design Criteria for Structures [19] and is utilized to compute the range of each variable.

Numerous numerical investigations have been carried out employing various constitutive models to analyze the behavior of piled raft systems, as reported by Sales et al. [20]. Additionally, Galliková and Rehman [21] conducted a settlement analysis of Neogene clay utilizing a modified hypoplastic-model-based finite element modeling. The hypoplastic model can provide greater accuracy than a simple Mohr–Coulomb or Cam-clay model for clayey soils. However, this study assumes a general sandy soil condition. In the numerical analysis, the applied ground is assumed to be sandy soil, and the Mohr–Coulomb model is employed, with the pile and raft modeled using the linear elastic model. The constitutive model and material properties utilized in the numerical analysis are outlined in

Table 2. Constitutive model and material properties for numerical analysis.

Material	Acronym	Unit	Value
			Loose Sand	Medium Sand	Dense Sand	Pile and Raft
Unit weight	γ	kN/m3	15	16	17
Young’s modulus	E	kPa	20 × 10−3	40 × 10−3	50 × 10−3	3.92 × 10−7
Poisson’s ratio	υ	-	0.3	0.3	0.3	0.15
Cohesion	c	kPa	0	0	0	-
Shear resistance angle	φ	° (deg)	35	40	45	-
Dilatancy angle	ψ	° (deg)	5	10	15	-
Rinter	-	-	0.8	0.6	0.5	-

.

Bolton [22] proposed the dilatancy angle as Equation (1), which is applied in this study, as shown in Table 2.

ψ = φ − 30°

(1)

Figure 1 shows a conceptual view of the numerical analysis, and Figure 2 summarizes the mesh generation of the numerical analysis.

The numerical modeling approach employed in this study encompasses a mesh consisting of over 68,500 individual elements. The model is subject to boundary conditions, where the x_min, x_max, y_min, and y_max values are fixed, while full fixity is imposed at z_min, and z_max is left unconstrained. Furthermore, the excavation of the tunnel is assumed to produce a −2% volume loss, while the groundwater conditions are disregarded in the model. To safeguard against any potential displacement of the piled raft affecting the boundary, the ground mass was secured at a distance of no less than 5B_R the width of the foundation.

2.2. Result of Numerical Analysis

This chapter presents an in-depth analysis of the pile axial stress distribution when tunnel excavation is carried out directly beneath the pile (OH = 0.0 D, OV = 0.5 D). The results are compared with the findings of Zheng et al. [11] based on field measurements and numerical analysis to verify the accuracy of the current study. Figure 3 illustrates the axial stress distribution of a pile at the pre-loading stage of tunnel excavation and at the end of tunnel excavation (TF = 0) in dense soil.

As demonstrated in Figure 3, the axial stress experienced by the pile decreases uniformly along its entire length during tunnel excavation from the pre-excavation stage in dense sand. The degree of reduction in the axial stress is found to be more pronounced in the central pile, as evident from Figure 3a, as compared with the left and right piles in Figure 3b,c. Furthermore, as depicted in Figure 4, the reduction in the axial stress is observed along the entire length of the pile in loose sand. However, in contrast with dense sand, the reduction in the axial stress is significantly greater in the middle and toe sections of the pile than in the upper portion. Additionally, while the central pile exhibits a greater reduction in axial force than the left and right piles in dense sand, the reduction in the axial force is observed across all piles in loose sand.

The changes in the pile axial stress are presented in Figure 5 and Figure 6, depicting the results for dense and loose sand, respectively. As demonstrated in Figure 5, the rate of change in the pile axial stress varies depending on the length of the pile in the center pile. Nevertheless, the axial stress in the center pile and at the top of the pile decreases overall during tunnel excavation. When the length of the pile is shorter, the middle and toe of the pile exhibit a more significant decrease compared with the top. However, for longer piles, the reduction in stress is less pronounced in the middle and toe. Notably, Figure 5b,c indicate that the axial stress increases at P_L = 10 and 20, with the decrease being insignificant at P_L = 35. Conversely, in loose sand, the decrease in pile axial stress is greater for shorter pile lengths and is greatest in the middle of the pile, as depicted in Figure 6. In addition, the left and right piles, shown in Figure 6b,c, respectively, exhibit a significant increase in pile axial stress at the top, followed by a rapid decrease. At PL = 35, however, there was no significant difference between the pile locations, and the overall rate of decrease was similar.

In the scenario where the tunnel excavation was carried out directly underneath the pile, the maximum pile axial stress was observed in the center pile, while the left and right piles exhibited a slight decrease or increase, in agreement with the findings of Zheng et al. [11]. Moreover, it can be inferred that the highest pile axial stress occurs at the top during reloading, implying that the results of the numerical analysis conducted in this study conform to the prevailing literature.

3. Prediction Model of Pile Axial Stress

3.1. Data Exploration and Preprocessing

This study utilizes a variety of variables and algorithms to derive and optimize a model capable of predicting the axial stress of a pile during tunnel excavation adjacent to a piled raft. The numerical analysis is used to identify the variables, and the resulting settlement is used as input data for machine learning. However, the allowable load applied to the piled raft is excluded from the input data for machine learning because it is expected to have a significant impact on the prediction model and, therefore, to predict the axial stress change accordingly.

In this study, the authors analyze the performance of the prediction models derived by varying the features of the input data.

Table 3. Features for machine learning in each case.

No.	Case 1	Case 2	Case 3	Case 4
1	EG	EG	EG	EG
2	OH	OH	OH	OH
3	OV	OV	OV	OV
4	PD	PD	PL/PD	PL/PD
5	PL	PL	APILE/ARAFT	APILE/ARAFT
6	BR	BR	PCTC	PCTC
7	PN	PN	TF	TF
8	PCTC	PCTC	PP	PP
9	TF	TF	PZ	PZ
10	PP	PP	Settlement
11	PZ	PZ
12	Settlement

shows that four cases are examined. In case 1, settlement is included as an independent feature among the variables applied to the numerical analysis, and in case 2, settlement is excluded from case 1 and applied to machine learning. In cases 3 and 4, the authors create new features by combining the variables considered in the numerical analysis. Specifically, A_PILE/A_RAFT is the cross-sectional area of the pile relative to the area of the raft in the piled raft, a variable created by combining P_D, P_N, and B_R from Table 1. Additionally, P_L/P_D is the length-diameter ratio of the pile, while P_P and P_Z indicate the location at which the pile axial stress is measured, with P_P divided into the left, center, and right based on the raft center and determined by applying the P_CTC value. P_Z is the depth at which the pile axial stress is measured. In case 3, the authors consider settlement as a feature for machine learning, whereas in case 4, settlement is excluded.

Table 3 provides a summary of the features considered for machine learning in each case.

Figure 7 shows that the null of the input data is missing in each case.

To identify the presence of any outliers within the data, the distribution of the dependent variable was subjected to examination, as illustrated in Figure 8. The analysis confirmed that no outliers are present within the data. While it is common for data obtained from field measurements to contain outliers, the present study did not extensively investigate the handling of outliers and instead relied on the results derived from numerical analysis.

The results of the correlation analysis between the features in each case are presented in Figure 9. It is important to note that correlation does not necessarily imply causality between the features, but rather a statistical relationship. A higher absolute value of the correlation coefficient indicates a stronger correlation, and the sign indicates the direction and proportionality of the relationship. Figure 9a illustrates the correlation analysis for case 1, where P_Z shows the highest correlation coefficient (0.49) with the target feature. Notably, as the value of P_Z decreases, the pile axial stress increases. Furthermore, E_G, P_L, and settlement exhibit correlation coefficients of 0.4, 0.36, and 0.26, respectively, with stress increasing as E_G and P_L increase and decreasing as settlement increases. Figure 9b depicts the correlation analysis for case 2, which is not significantly different from the analysis in case 1 except for the absence of the settlement feature.

3.2. Selection of Algorithms

In this study, the authors conduct a comparative analysis of the predictive models of pile axial stress derived in each case in terms of their performance, utilizing various machine-learning algorithms. Specifically, eight algorithms are considered, namely, Linear Regression (LR), Lasso (L1), Ridge (L2), Multi-layered Perceptron (MLP), Support Vector Machine (SVR), Random Forest (RF), Extreme Gradient Boosting (XGB), and Light Gradient Boosting (LightGBM).

LR is a straightforward, rapid, and comprehensible algorithm used for the linear prediction of regression problems. However, it has the drawback of being less accurate. L1 and L2 are algorithms designed to address the shortcomings of LR by imposing additional constraints to prevent overfitting of the predictive models derived through LR. L1 and L2 are restricted via L1-norm and L2-norm penalties, respectively. A detailed explanation of the LR, L1, and L2 algorithms can be found in the following references: [23], [24], and [25].

MLP is a multilayer neural network structure with multiple perceptrons in various layers, and is a fundamental algorithm of artificial neural networks, containing one or more hidden layers between the input and output layers. For further details on the MLP algorithm, the reader is referred to [26].

SVR is one of the most popular algorithms used in machine learning and separates data by choosing the hyperplane that is furthest from the data. It is less sensitive to erroneous data, less prone to overfitting, and more user-friendly than neural network algorithms such as MLP. However, it has the drawbacks of taking longer to learn and being less interpretable. A detailed explanation of the SVR algorithm can be found in [27].

Random Forest is one of the representative ensemble algorithms that combines multiple decision-tree algorithms to derive a predictive model. It uses a bagging method for training and has the advantages of good generalization performance, ease of parameter adjustment, and poor overfitting. However, it also has the disadvantages of difficult analysis for individual trees and poor model performance when the dimensionality is large and the amount of data is not significant. A description of RF can be found in [28].

XGB is a boosting-based algorithm that has won several competitions and is one of the ensemble algorithms that, like RF, uses decision trees. For a description of XGB, the reader should refer to [29].

LightGBM, like XGB and RF, is an ensemble algorithm that utilizes a decision tree. Unlike XGB, which expands the tree in a level-wise manner, LightGBM expands the tree leaf-wise. This reduces the number of losses when expanding the same leaf, and the algorithm has the advantage of being faster than XGB. For an explanation of LightGBM, the reader can refer to [30].

In machine learning, the data are typically divided into training, validation, and test data in a ratio of 5:3:2 to prevent overfitting, with only the training data being used for algorithm selection. The learning results of each algorithm were evaluated using the root-mean-square error (RMSE), and the algorithm results are summarized in

Table 4. RMSE of the training of each algorithm.

	Case 1	Case 2	Case 3	Case 4
LR	763.97	769.57	744.39	754.37
L1	767.49	772.81	744.37	754.34
L2	769.18	774.53	744.61	754.87
MLP	997.53	1100.50	1074.01	994.62
SVR	1352.94	1352.64	1352.45	1352.22
RF	207.17	289.75	204.55	291.34
XGB	187.29	276.95	185.06	275.98
LightGBM	184.28	244.13	186.97	241.69

.

The results reveal that the algorithms RF, XGB, and LightGBM exhibited the best performance, while MLP and SVR, which are based on artificial neural networks (ANN), demonstrated a lower performance. Therefore, hyper-parameter tuning was performed only on RF, XGB, and LightGBM to enhance the predictive model performance.

3.3. Hyper-Parameter Tuning

Hyper-parameters constitute the intrinsic parameters of an algorithm that can be fine-tuned to optimize the performance of a prediction model. As each algorithm has its unique set of hyper-parameters, modifying them inevitably impacts the predictive model’s performance. The hyper-parameters for the aforementioned algorithms are listed in

Table 5. Hyper-parameters of each algorithm.

	RF		XGB		LightGBM
n_estimators	426	n_estimators	96	n_estimators	140
max_depth	7	max_depth	8	max_depth	14
min_samples_split	10	learning_rate	0.15	learning_rate	0.06
max_leaf_nodes	20	colsample_bytree	0.93	colsample_bytree	0.95
random_state	21	subsample	0.78	subsample	0.07
		reg_alpha	41.21	reg_alpha	47.2
		reg_lambda	35.17	reg_lambda	11.06
		gamma	0.05	random_state	23
		random_state	16

for reference.

A description of the hyper-parameters considered for RF, XGB, and LightGBM can be found in the following references: [28], [29], and [30]. Hyper-parameter tuning involves optimizing the parameters to achieve the best possible performance of the predictive model applied with the algorithm. In this study, the Optuna framework was employed to search for the optimal value by adjusting the hyper-parameters in each trial within the specified range or set list. Hyper-parameter tuning was carried out using the validation data from the previously divided data in the algorithm selection phase. The results of the hyper-parameter tuning and the resulting predictive models for each case are summarized in

Table 6. Results of hyper-parameter tuning for each case.

	Case 1	Case 2	Case 3	Case 4
RF	328.80	331.02	301.02	313.47
XGB	176.87	233.88	178.26	234.45
LightGBM	173.61	232.25	173.34	232.89

.

After tuning the hyper-parameters using the validation data, LightGBM showed the best performance.

3.4. Performance of Prediction

As previously discussed, the test data were predicted using the prediction model obtained from the algorithm selection process using the training data and the optimization of the hyper-parameters through the validation data. It is crucial to note that the test data must not be exposed to the algorithm during the training and validation phases to prevent overfitting. The results of the hyper-parameterization, as elaborated in Section 3.3, are succinctly presented in Figure 10, Figure 11, Figure 12 and Figure 13.

The results of the comparison of the models in terms of performance, as predicted for each case and algorithm using the RMSE, are summarized in

Table 7. Performance of prediction model.

	Case 1	Case 2	Case 3	Case 4
RF	332.51	381.24	307.61	327.62
XGB	203.08	247.21	199.63	265.41
LightGBM	200.82	262.23	197.65	257.51

.

Table 7 presents the performance results of LightGBM and XGB, both of which displayed satisfactory training and validation performance, with case 3 showing the most accurate predictions. Comparing cases 1 and 3 with cases 2 and 4, it is evident that adjusting the features in cases 3 and 4 resulted in better performance than applying the same features as in the numerical analysis, as observed in cases 1 and 2. Although the final prediction is slightly higher than the RMSE of the model obtained from the training and validation, the difference is negligible, and there is no loss of confidence in the model due to overfitting. When calculating the error of the best-performing LightGBM prediction model in case 3 with respect to the axial force of the pile, the error ranges from approximately 39 kN to 155 kN for pile diameters between 0.5 m and 1.0 m.

4. Explanation of the Prediction Models Using Explainable AI

4.1. Shapley Value

One commonly used technique to analyze the relationship between a predictive model and its variables is feature importance. However, this method can indicate only the degree to which each independent feature influences the predicted value, without revealing how much the predicted value changes as each feature is varied. Furthermore, feature importance is applicable only to tree- or boosting-based ensemble algorithms, making it challenging to discern the relationship between individual features and predicted values in models derived through other algorithms, particularly those based on artificial neural networks.

To address this limitation, Shapley values can be employed. These values rely on game-theory principles to calculate the contribution of each player to the gain, and are determined by constructing a combination of attributes to assess the importance of a single attribute and averaging the change in its presence or absence. In recent years, the SHAP framework has leveraged Shapley values to enhance the explanatory power of predictive models in a variety of fields related to machine learning and deep learning [31].

4.2. Analysis of the Prediction Models

In this chapter, the authors employ the feature importance and SHAP summary plot results to illustrate the predictive models derived in each case. Figure 14 displays the importance of the features in each case, quantifying the impact of each feature on the predictive model. Specifically, Figure 14a showcases the impacts of independent features on the predictive models derived by the RF, XGB, and LightGBM algorithms in case 1. For the best-performing LightGBM, the settlement exhibits the highest impact at 43.17%, followed by P_Z and P_L at 9.94%, T_F at 6.44%, O_H at 6.00%, O_V at 5.92%, P_P at 5.40%, P_CTC at 5.24%, EG at 3.44%, P_D at 2.55%, and P_N at 1.95%; B_R does not affect the prediction model at all. The feature importance in case 2 is illustrated in Figure 14b, where P_L has the highest impact at 22.02%, followed by P_CTC at 15.26%, P_P at 11.42%, E_G at 10.09%, P_D at 10.04%, P_Z at 8.27%, T_F at 6.62%, P_N at 6.55%, and O_V and O_H at 4.40% and 3.99%, respectively; B_R has a negligible impact at 1.35%. Case 3 is summarized in Figure 14c, where the settlement has the highest impact at 42.01%, followed by P_Z at 10.92%, P_L/P_D at 8.87%, T_F at 6.99%, P_P at 5.74%, O_V at 5.60%, P_CTC at 5.92%, O_H at 5.55%, A_PILE/A_RAFT at 5.09%, and EG at 3.30%. Figure 14d shows the feature importance in case 4, where P_CTC has the highest impact at 19.02%, followed by A_PILE/A_RAFT at 16.03%, P_P at 15.53%, P_Z at 14.54%, T_F at 9.29%, E_G at 8.07%, P_L/P_D at 7.96%, and O_H and O_V at 4.81% and 4.75%, respectively.

The analysis of the feature importance reveals that settlement is of significant importance in cases 1 and 3, where it is included as an independent variable. Meanwhile, P_CTC has the highest feature importance in cases 2 and 4, where settlement is not considered. However, while the feature importance results can quantify the impacts of features on the predictive value of a target feature, they do not provide insight into whether the prediction value increases or decreases as the number of features increases. Therefore, the impact of each feature on the prediction model, as measured by the SHAP value, is analyzed in the summary plot presented in Figure 15.

In the summary plot, a large absolute value of the SHAP value indicates a large contribution to the predicted value. A positive SHAP value indicates that as the feature value increases, the prediction value also increases, and as the feature value decreases, the prediction value decreases.

Figure 15a illustrates the contributions of the features in the model derived by the best-performing LightGBM algorithm in case 1 as a function of the SHAP value, with the predicted values shown in order of magnitude. The feature with the largest contribution is P_L, with the SHAP value increasing as the feature value increases. This trend indicates that the pile axial stress increases with the length of the pile, likely due to the larger allowable load of long piles compared with short piles. Next is P_Z, which demonstrates that the pile axial stress increases as the feature value decreases, meaning that the pile axial stress is greater at the top of the pile. E_G indicates that a larger relative density of the ground results in a greater allowable load and, therefore, greater axial stress on the pile. P_D shows that the pile axial stress decreases as the feature value increases, indicating a small change in the pile axial stress due to the tunneling. Settlement is found to decrease as the feature value increases, as the settlement caused by the tunnel excavation reduces the pile axial stress. P_P, which represents the position of the pile, shows that the smaller the feature value, the larger the SHAP value, meaning that the greater distance from the tunnel side results in greater pile axial stress. P_CTC shows that the smaller the feature value, the smaller the SHAP value, suggesting that the closer the piles are in terms of spacing, the smaller the pile axial stress, due to the cluster pile effect. The feature P_N shows that the larger the feature value, the smaller the SHAP value, indicating that a larger number of piles results in smaller pile axial stress, as more piles share the allowable load. However, T_F, B_R, O_H, and O_V do not produce meaningful SHAP values in terms of their impacts on the predictive model.

Figure 15b shows the SHAP values of the prediction model derived by LightGBM in case 2, and the analysis reveals a similar trend to that observed in case 1. However, the B_R in case 2, which was not analyzed in case 1, exhibits a more meaningful SHAP value than in case 1, with the SHAP value decreasing as the feature value increases. This trend is due to the fact that the larger the foundation width, the more of the allowable load is supported by the raft, resulting in less pile axial stress.

The SHAP values in case 3, summarized in Figure 15c, reveal that the smaller the feature value of P_L/P_D, the smaller the SHAP value, which is consistent with cases 1 and 2. Additionally, the larger the A_PILE/A_RAFT, the smaller the SHAP value, with a larger A_PILE indicating a smaller A_RAFT. A_PILE is a function of P_D and P_N, while A_RAFT is a function of B_R, and larger P_D and P_N values result in smaller pile axial stress.

5. Discussion

This study employs numerical analysis and machine-learning techniques to develop a model for predicting pile axial stress during tunnel excavation adjacent to a piled raft. The numerical analysis considers various design parameters for the piled raft and tunneling, and uses settlement and pile axial stress results from multiple cases as input data for machine learning. The impact of each independent variable on the prediction model is analyzed using feature importance and SHAP values.

During the training phase, RF, XGB, and LightGBM algorithms were applied, with LightGBM exhibiting the best performance. This trend is observed consistently across all cases, indicating that these three algorithms are the most reliable for deriving a predictive model of pile axial stress. In the validation phase, LightGBM proved to be the most accurate and exhibited the lowest RMSE in all cases, as compared with RF and XGB. In the prediction phase, LightGBM performed the best in all cases, followed by cases 3, 1, 4, and 2. For cases 1 and 2, the same variables used in the numerical analysis were employed, while a combination of features was used for cases 3 and 4. The model that does not account for settlement performed the best, but overall, the model that incorporates settlement as an independent feature performed well.

This finding is supported by the feature importance and SHAP value analyses, which indicate that settlement is a significant variable to consider. According to the SHAP value analysis, factors that reduce the pile axial stress due to tunnel excavation include pile top, long pile, loose ground, large pile diameter, large foundation settlement, narrow pile spacing, close pile–tunnel distance, large number of piles, and large foundation width.

6. Conclusions

The present study reveals that while developing a prediction model using machine-learning or deep-learning techniques, it is crucial to take into account not only the algorithm but also the method used for processing the input data, especially when a non-numerical laboratory model or field monitoring data are used. Additionally, several relevant variables must be considered to predict the pile axial stress during the excavation of a tunnel adjacent to a piled raft in practical applications. To obtain a more precise estimation, it is recommended to incorporate measured settlement data of the foundation into the analysis. Consequently, the authors intend to develop a prediction model by processing diverse types of data and enhance its reliability by testing it with laboratory model and field monitoring data.