Explainable GAN-Augmented MLP for Soil Resilient Modulus Prediction

1  Introduction

The resilient modulus (MR) of subgrade soils is a critical mechanical parameter for the design and performance evaluation of pavement structures. Conventionally, MR is determined through repeated load triaxial (RLT) tests in laboratory settings [1]. However, this approach is often criticized for being time-consuming and economically inefficient [2,3]. As a result, numerous researchers have focused on developing alternative models to estimate MR indirectly, thereby improving the efficiency of pavement design processes [3].

Various soil properties, including water content, unit weight, plasticity index, and liquid limit, have been shown to significantly influence the value of MR [4,5]. Building on these relationships, several empirical and semi-empirical models have been proposed. For instance, Mehrotra et al. developed a predictive method using basic soil properties while accounting for moisture and stress conditions [6]. Similarly, Sadrossadat et al. developed an indirect estimation method considering soil properties and stress state, showing superior predictive ability compared to traditional models [7]. However, laboratory-derived parameters may not fully capture field conditions, highlighting the importance of in situ tests for more reliable MR predictions [8,9]. For example, Mousavi et al. utilized dynamic cone penetrometer (DCP) data to characterize in situ soil properties and predict MR across varying stress states [10].

Despite these advancements, empirical and semi-empirical methods still rely heavily on predefined functional forms, manual parameter calibration, and domain-specific assumptions, which may limit their adaptability when soil conditions, stress states, or testing environments vary. To overcome these limitations, machine learning (ML) methods have been increasingly adopted for MR prediction because of their ability to capture nonlinear relationships between soil properties and mechanical responses without explicitly prescribing the model form. Recent studies have applied various ML algorithms to MR estimation. For example, Pahno et al. [11] compared regression tree, Random Forest (RF), extreme gradient boosting (XGBoost), and multiple linear regression models using routinely collected soil properties, and reported that tree-based ensemble models achieved better predictive performance than conventional regression models. Kardani et al. [12] developed several standalone and ensemble ML models based on a large database of compacted subgrade soils, further demonstrating the potential of ensemble learning for MR prediction. Khawaja et al. [13] compared genetic programming, multi-expression programming, and artificial neural network (ANN)-based models for forecasting subgrade soil MR. More recently, Bardhan et al. [14] improved the predictive capability of RF by coupling it with the grey wolf optimizer and compared the optimized model with several regression and neural-network-based models. These studies indicate that ML methods provide flexible and effective alternatives to empirical equations for modeling the complex relationship between soil properties and MR.

Beyond conventional ML methods, deep learning (DL) has achieved remarkable success in many fields because of its strong capability to learn complex nonlinear representations from data. In geotechnical engineering, DL models have also been increasingly adopted for solving prediction and classification problems involving complex soil behavior [15–17]. For MR prediction, neural-network-based models have been successfully applied to both cohesive and non-cohesive soils, achieving promising predictive accuracy [18]. Compared with conventional ML models, DL models offer greater flexibility in representing interactions among multiple soil variables and mechanical responses. A concise comparison of representative data-driven studies on MR prediction is provided in Table 1.

As summarized in Table 1, existing data-driven studies have advanced MR prediction through model selection, ensemble learning, and metaheuristic optimization. However, these studies have mainly focused on improving the predictive model itself, while the limited availability of site-specific geotechnical data remains a fundamental challenge. This issue is particularly important for DL models, which have shown great potential in MR prediction [18] but generally require sufficient and representative training samples to learn complex nonlinear relationships. When the available dataset is small, DL models may suffer from overfitting and may not always outperform conventional ML methods [13]. In addition, their predictive performance is highly dependent on hyperparameter selection, and manual tuning is usually labor-intensive and suboptimal [19]. Furthermore, DL models are often regarded as “black-box” models, which limits their interpretability and practical acceptance in engineering applications [20,21]. Therefore, improving small-sample MR prediction requires not only an effective predictive model, but also a framework that can address data scarcity, hyperparameter optimization, and model interpretability simultaneously.

Data augmentation provides a potential strategy for alleviating data scarcity in small-sample modeling. Among different data augmentation methods, Generative Adversarial Networks (GANs) have demonstrated considerable potential in synthesizing realistic data samples that mimic real-world distributions. A GAN operates through a competitive learning process between two neural networks, namely a generator and a discriminator [22]. This capability enables GANs to augment datasets by generating additional synthetic samples, thereby improving the training of data-driven models under limited-data conditions. GANs have achieved significant success in areas such as image generation [23] and data enhancement [24]. In geotechnical engineering, data augmentation strategies have also been explored to improve model performance under limited-data conditions. For example, Jiang et al. [25] developed a data-augmentation-assisted CNN framework for probabilistic slope stability analysis in spatially variable soils, while Chen et al. [26] used Wasserstein GANs to enhance small-sample soil liquefaction datasets. However, the application of GAN-based data augmentation to soil property modeling, particularly for MR prediction, remains underexplored. This highlights the need for dedicated research on GAN-based augmentation strategies tailored to small-sample soil datasets.

To address these challenges, this study develops a GAN-Optuna-multilayer perceptron (MLP) framework for MR prediction under small-sample conditions. The framework uses GAN-based data augmentation to enrich the limited training dataset, Optuna-based optimization to determine suitable MLP hyperparameters, and Shapley Additive Explanations (SHAP) analysis to explain the influence of individual input variables on the predicted MR. This integrated strategy is intended to improve the predictive performance and interpretability of data-driven MR modeling.

2  Materials and Methods

This section describes the materials and methods used to develop the proposed GAN-Optuna-MLP framework for predicting soil MR. First, the dataset adopted in this study is introduced, including the soil type, data sources, input variables, and target variable. Then, the main methodological components, including GAN, MLP, Optuna, SHAP, k-fold cross-validation, and evaluation metrics, are presented. Finally, the implementation procedure of the proposed framework is described.

2.1 Dataset

DL-based algorithms are gaining prominence in modeling the MR of soils. For these algorithms, the selection of appropriate training data is crucial. The database used in this study was derived from the clay database reported by Liu et al. [3], which focuses on cohesive subgrade soils. The soil samples were collected from 16 sites in Jiangsu Province, China, and the subsoil strata mainly consist of Quaternary soft to stiff clayey soils and silty clay soils. Therefore, the proposed model was developed and evaluated for cohesive subgrade soils within the range of the adopted database. The database includes laboratory test results and cone penetration test results, which provide the MR value, cone tip resistance (qc), sleeve friction resistance (fs), water content (w), and dry density (γd) of the soil. The database is relatively small, containing only 124 samples. This limited sample size provides a suitable case for evaluating the effectiveness of the proposed GAN-based data augmentation method.

2.2 GANs

GANs are a class of generative learning algorithms inspired by two-player zero-sum games from game theory. A GAN consists of two main components: a generator and a discriminator. The generator’s role is to create synthetic data, while the discriminator distinguishes between synthetic and real data. These two components are trained concurrently in a competitive manner, where the generator aims to produce data that is indistinguishable from real data, gradually improving over time.

The training of a GAN involves an alternating optimization process, which is divided into two stages. In the first stage, the discriminator is fixed, and the generator’s parameters are adjusted to generate synthetic samples that are difficult for the discriminator to distinguish. In the second stage, the generator is fixed, and the discriminator’s parameters are optimized to improve its ability to distinguish real samples from generated samples. The generator’s objective is to minimize the loss function described in Eq. (1):

VG=Ez∼Pz(z)[log⁡{1−D(G(z))}](1)

Conversely, the discriminator aims to maximize the loss function given by Eq. (2):

VD=Ex∼Pdata(x)[log⁡D(x)]+Ez∼Pz(z)[log⁡{1−D(G(z))}](2)

By combining Eqs. (1) and (2), the objective function for the GAN training process can be derived as:

Ex∼Pdata(x)[log⁡D(x)]+Ez∼Pz(z)[log⁡{1−D(G(z))}](3)

where E denotes the mathematical expectation; x represents a real sample from the training dataset; z is a latent noise vector sampled from the prior distribution Pz(z); Pdata(x) denotes the distribution of the real data; G represents the generator; G(z) is the synthetic sample generated by the generator; D represents the discriminator; D(x) is the discriminator output for a real sample; and D(G(z)) is the discriminator output for a generated sample, indicating the probability that the generated sample is regarded as a real sample by the discriminator.

2.3 MLP

An MLP is a traditional DL architecture that forms a network of artificial neurons arranged in a specific topology. This architecture includes input, hidden, and output layers. The input layer receives raw data, hidden layers perform complex data transformations and feature extraction, and the output layer delivers the final computational results. Information flows through the network via weighted connections between layers, where each connection’s weight parameter indicates the strength of information transfer. Additionally, each neuron employs an activation function to introduce nonlinearity, enabling the network to learn and model complex functional relationships.

2.4 Optuna

Optuna is an automatic hyperparameter optimization framework proposed by Akiba et al. [27]. It adopts a define-by-run design, allowing users to dynamically construct the hyperparameter search space and optimize a user-defined objective function. This feature makes Optuna suitable for DL models, whose performance is strongly affected by architecture- and training-related hyperparameters, such as the learning rate, number of hidden neurons, and regularization parameters.

In Optuna, each hyperparameter combination is evaluated through a trial. During each trial, a candidate set of hyperparameters is sampled, the model is trained using this configuration, and the validation loss or another predefined metric is returned as the objective value. Based on the results of previous trials, Optuna adaptively suggests new candidate hyperparameters and searches for the configuration that minimizes or maximizes the predefined objective function. The tree-structured Parzen estimator (TPE), which is commonly used in Optuna, models the distributions of hyperparameter configurations associated with better and worse objective values separately, and then samples new candidates from regions that are more likely to improve the objective value [27].

Optuna-based hyperparameter optimization has been successfully applied in various data-driven prediction tasks, including pavement material characterization, tunnel-induced settlement prediction, and other geotechnical engineering problems [28,29]. These applications indicate that Optuna provides an efficient and reliable alternative to manual trial-and-error tuning, particularly for ML models with multiple interacting hyperparameters.

In this study, Optuna was used to optimize the hyperparameters of the prediction models. For the MLP model, Optuna evaluated different hyperparameter combinations through multiple trials, and the combination with the lowest validation loss was selected for final model training. To provide a fair comparison, RF and XGBoost were also introduced as benchmark models and optimized using the same Optuna framework under the same training/testing split.

2.5 SHAP

DL models, currently the most widely used nonlinear models, are applied across diverse domains. However, they are often considered “black-box” models due to their limited interpretability. SHAP [30], based on game theory, provides an interpretable method to explain the output of any ML model [31]. The SHAP method can measure the impact of features on the dependent variable in ML models, assigning an importance value (SHAP value) to each feature. Its expression is:

φi=∑S⊆N∖{i}|S|!(|N|−|S|−1)!|N|

where: N is the set of all features. S is a subset of N that does not include feature i. f(S) is the model prediction when only the feature set S is involved. f(S∪{i}) is the model prediction when feature i is added to the feature set S. The combinatorial coefficient |S|!(|N|−|S|−1)!|N|! is used to weight the influence of each subset of features.

2.6 K-Fold Cross-Validation

In model training, particularly when data availability is limited, cross-validation is the predominant method for determining model hyperparameters. The most commonly used technique is k-fold cross-validation. This method involves randomly dividing the original training dataset into k subsets. In each iteration, k − 1 subsets are used for training, and the remaining subset is used as the validation set. This process is repeated k times, with the model’s performance evaluated based on the average prediction error across all subsets. In this study, we implemented 5-fold cross-validation to enhance the robustness of the model, using the average loss across the validation sets (Eq. (5)) as the objective for Optuna-based hyperparameter optimization.

Lavg=15∑i=15Li(5)

where Li represents the loss for the i-th fold, and Lavg denotes the final averaged loss.

2.7 Evaluation Metrics

In this study, the coefficient of determination (R2), mean absolute error (MAE), and mean absolute percentage error (MAPE) were used to assess the performance of different models. The definitions of R2, MAE, and MAPE are as follows:

R2=1−∑i=1N(yip−yi)2∑i=1N(yi−y¯)2(6)

MAE=1N∑i=1N|yip−yi|(7)

MAPE=1N∑i=1N|yip−yi|yi×100%(8)

where N represents the number of samples, yi represents the target data of the i-th set of samples, yip represents the predicted value of the i-th sample, and y¯ represents the mean value of all measured samples.

2.8 GAN-Optuna-MLP

Fig. 1 illustrates the framework of the proposed GAN-Optuna-MLP model. Initially, the data is partitioned into training and testing sets. The training set is processed by a GAN model to generate an augmented dataset, which is then split into new training and validation sets using k-fold cross-validation. The model hyperparameters were optimized using the Optuna framework based on the validation performance obtained from k-fold cross-validation. The optimal model and its hyperparameters are determined through the Optuna optimization process. These hyperparameters are subsequently employed to train DL models, which are then used to predict outcomes.

Figure 1: GAN-Optuna-MLP based MR prediction framework.

It should be noted that the GAN was not trained on the full dataset. In this study, the GAN was trained only on the training set after the initial train–test split. The generated synthetic samples were combined with the original training samples to construct the augmented training dataset. The independent testing set remained unseen during GAN training, data augmentation, and hyperparameter optimization; therefore, data leakage was avoided. The held-out testing set was used to simulate unseen samples for evaluating the generalization ability of the proposed model.

2.9 Training of the GAN-Optuna-MLP Model

This study utilized a GAN model for data generation, comprising a generator and a discriminator. Table 2 outlines the generator’s architecture with three layers: an input layer corresponding to a 32-dimensional latent noise vector, a hidden layer mapping 32 units to 64 units via the ReLU activation function, and an output layer of 5 units employing the Sigmoid activation function. Table 3 presents the discriminator’s architecture, which includes a 5-unit input layer, a 64-unit hidden layer using ReLU, and a 1-unit output layer with Sigmoid activation.

The GAN-generated synthetic samples were combined with the original training samples and then used to train the GAN-Optuna-MLP model. In DL models, selecting the correct hyperparameters is essential for achieving reliable results. These hyperparameters, optimized via Optuna, are outlined in Table 4. The model employs the Adam optimizer to update weights, with a batch size of 32 and a learning rate of 0.00188. Training the model for 100 epochs has been found to yield optimal results. Additionally, weight decay is utilized to prevent overfitting.

To explicitly evaluate the model performance during hyperparameter optimization, the validation results obtained from 5-fold cross-validation were further reported. In addition to Optuna-MLP and GAN-Optuna-MLP, RF and XGBoost were introduced as benchmark models. Their hyperparameters were also optimized using Optuna under the same training/testing split to ensure a fair comparison. Therefore, these benchmark models are denoted as Optuna-RF and Optuna-XGBoost in Table 5.

During the Optuna optimization process, four folds were used for model training, while the remaining fold was used for validation. The validation performance was calculated as the mean and standard deviation over the five folds, as shown in Table 5. This validation procedure was used for hyperparameter selection, while the independent testing set was kept unseen during model training, GAN-based data generation, and hyperparameter optimization.

As shown in Table 5, the proposed GAN-Optuna-MLP model achieved the best validation performance among all models. Compared with Optuna-MLP, the GAN-Optuna-MLP model showed higher validation accuracy and lower prediction errors, indicating that GAN-based data augmentation effectively improved the predictive performance of the MLP model under small-sample conditions. In addition, the relatively small standard deviations suggest that the proposed model maintained stable performance across different validation folds.

Fig. 2 displays the loss trajectories for both training and testing data across epochs. The accuracy of the model can be inferred from the reduction in training loss over time. As depicted, the loss decreases progressively and stabilizes with increasing epochs, indicating effective learning. By the end of the training period, the model demonstrates high accuracy on both training and testing data.

Figure 2: Evolution of loss value during the training process.

3  Results and Discussion

This section presents the results and discussion of the proposed framework. First, the quality of the GAN-generated data is evaluated using distributional statistics and correlation analysis. Then, the predictive performance of the proposed GAN-Optuna-MLP model is compared with those of Optuna-MLP, Optuna-RF, and Optuna-XGBoost benchmark models to assess the effectiveness of GAN-based data augmentation under small-sample conditions. A sensitivity analysis is subsequently conducted to examine the influence of GAN randomness on model performance. Finally, SHAP analysis is used to interpret the contribution of each input variable to MR prediction.

3.1 Evaluation of GAN-Generated Data

This study utilizes a GAN to generate synthetic data. A total of 200 synthetic samples were generated using the trained GAN model. Kernel density estimations (KDE) were employed to evaluate the similarity between the generated and actual data. Fig. 3 displays a KDE comparison of various features, illustrating that the synthetic data closely matches the distribution of the actual data, effectively demonstrating the GAN’s capability to replicate the original dataset’s statistical properties. Generally, the GAN-generated data aligns well with the KDE profiles of the actual data, indicating successful feature replication.

Figure 3: KDE comparison between GAN-generated and real data for different features: (A) MR; (B) qc; (C) fs; (D) w; (E) γd.

Table 6 presents the descriptive statistics for both the synthetic and real datasets. The results indicate that the distribution of the generated data is generally consistent with that of the actual data. Notably, relatively large discrepancies between the synthetic and real data are observed in the maximum and minimum values of certain features. This can be attributed to the limited sample size in the extreme ranges of the distribution, which potentially constrains the GAN model’s ability to accurately replicate the real data in these regions.

Fig. 4 illustrates the correlation matrices among MR, qc, fs, w, γd for the original and GAN-generated datasets. Since the objective of this study is to predict MR, the analysis mainly focuses on the correlations between MR and the input variables. In the original dataset, MR shows a strong positive correlation with qc (0.79), a moderate positive correlation with fs (0.47) and γd (0.46), and a clear negative correlation with w (−0.68). These trends are consistent with geotechnical understanding: higher cone resistance, sleeve friction, and dry density generally indicate stronger or denser soil conditions, whereas higher water content tends to reduce soil stiffness.

Figure 4: Correlation matrix of MR and input features: (A) original dataset; (B) GAN-generated dataset.

Similar correlation patterns are observed in the GAN-generated dataset. The correlation coefficients between MR, qc, fs, w, γd are 0.82, 0.50, −0.65, and 0.50, respectively. The differences between the original and generated datasets are relatively small, indicating that the GAN-generated data preserve the main correlation structure between MR and the influencing factors. The preservation of these physically consistent correlations suggests that the synthetic samples are not only statistically similar to the real data, but also physically plausible for cohesive subgrade soils. This further supports the reliability of the generated dataset for subsequent model training.

Taken together, the KDE curves, descriptive statistics, and correlation matrices indicate that the GAN-generated data generally preserve the main distributional characteristics and correlation structure of the original dataset. This is important for the proposed framework because the original MR database is limited in size, and insufficient samples may restrict the training of DL models and increase the risk of overfitting. By providing additional samples that follow similar statistical patterns to the real data, the synthetic dataset enriches the training information available for the MLP model and helps it learn the nonlinear relationships between MR and the input variables under small-sample conditions.

3.2 Validation and Comparison of Prediction Results

Fig. 5 displays the actual measured MR alongside the DL model predictions for both training and testing sets. The figure illustrates minimal errors between predicted outcomes and actual measurements, confirming the model’s high predictive accuracy.

Figure 5: Predicted and measured MR for all samples.

Fig. 6 compares the measured and predicted MR values of the GAN-Optuna-MLP and Optuna-MLP models for the training and testing sets. The values align closely with a slope of one, indicating robust model performance. Despite general alignment, some deviations occur, likely due to inherent data inaccuracies from field and lab tests, thus precluding perfect correspondence between predictions and actual results. Notably, the GAN-Optuna-MLP model’s predictions align most closely with the ideal line, suggesting superior accuracy among the evaluated models.

Figure 6: Comparison between measured and predicted MR values for GAN-Optuna-MLP and Optuna-MLP: (A) training set; (B) testing set.

To further evaluate the predictive performance of the proposed model, the optimized benchmark models described above were compared with Optuna-MLP and GAN-Optuna-MLP on the independent testing set. All models were evaluated using the same input variables, data partitioning strategy, and performance metrics.

As shown in Table 7, the Optuna-RF and Optuna-XGBoost models outperform the Optuna-MLP model without GAN-based data augmentation, indicating that tree-based ensemble learning models showed better generalization than the conventional MLP model on the original small-size dataset in this study. This result also suggests that the conventional MLP model may have limited generalization ability when trained only on the original small dataset.

After introducing GAN-based data augmentation, the predictive performance of the MLP model is substantially improved. Compared with Optuna-MLP, GAN-Optuna-MLP increases the testing R2 by 6.9%, reduces the testing MAE by 24.9%, and decreases the testing MAPE by 25.1%. These improvements demonstrate that the GAN-generated samples can enrich the limited training dataset and enhance the generalization ability of the MLP model under small-sample conditions.

Compared with the optimized ensemble learning models, the proposed GAN-Optuna-MLP model achieves the same testing R2 of 0.93, while yielding the lowest testing MAE of 2.72 MPa and MAPE of 5.70%. This indicates that, although RF and XGBoost are strong benchmark models for small-size tabular datasets, the proposed GAN-Optuna-MLP framework can achieve competitive or slightly better testing performance by alleviating the data scarcity problem of MLP-based prediction. Overall, these findings further support the effectiveness of GAN-based data augmentation for soil MR prediction.

3.3 Robustness Analysis under Random Initialization

To further evaluate the impact of randomness on prediction performance, a sensitivity analysis was performed through multiple independent training runs. In this analysis, the original training/test split was kept unchanged, while the GAN and the MLP model were reinitialized in each run. For each run, the GAN was retrained to generate an augmented dataset, and the resulting GAN-Optuna-MLP model was evaluated on the same independent test set.

The prediction performance was reported as the mean value and standard deviation over repeated runs, as shown in Table 8. The GAN-Optuna-MLP model achieved R2 = 0.93 ± 0.005, MAE = 2.73 ± 0.06 MPa, and MAPE = 5.67 ± 0.03%. The small standard deviations indicate that the proposed framework is not highly sensitive to stochastic training effects, including GAN-based data generation and model initialization, and can maintain stable prediction performance under different random initializations.

3.4 Interpreting the DL Model with SHAP Method

The DL model proposed in this study successfully predicted MR. However, due to the inherent characteristics of DL models, the exact reasoning behind the predictions remains unknown. This black-box nature of DL models poses a challenge, as understanding the contribution of each feature to the model’s predictions is crucial for validating the model and ensuring its reliability.

To address this issue, the SHAP method was employed to interpret the DL model’s predictions. The SHAP method allows us to obtain SHAP values, which reflect the differences between features and provide a unified measure of feature importance. This method calculates SHAP values for various features, allowing for an understanding of how the model’s predictions are influenced by each feature. The resulting SHAP values of the considered features are presented in Fig. 7.

Figure 7: SHAP values for different features.

The dominant role of qc is consistent with geotechnical understanding. Cone tip resistance is an in-situ mechanical index that reflects the overall resistance of soil during penetration and is closely associated with soil strength, stiffness, stress state, and stratigraphic variability. Therefore, qc can provide direct information on the mechanical response of subgrade soils under field conditions. This finding is consistent with Liu et al. [3], who reported that qc was the most informative single parameter for MR estimation among the considered penetration-related and physical indices.

The second most influential variable is water content w. For cohesive subgrade soils, water content strongly affects the resilient response by changing matric suction, effective stress, and the bonding/frictional resistance within the soil skeleton. An increase in water content generally weakens the soil structure and reduces stiffness, particularly for unsaturated fine-grained soils. Therefore, the relatively high SHAP importance of w is physically reasonable and indicates that the model captures the moisture sensitivity of MR.

The contribution of sleeve friction fs is lower than that of qc. Although fs reflects the frictional behavior along the cone sleeve, it is partly related to qc and may provide partially overlapping information. As a result, its independent contribution to the prediction of MR is smaller than that of qc. This interpretation is also consistent with the multivariate correlation analysis of Liu et al. [3], which suggested that prediction accuracy is not necessarily improved when strongly correlated input variables are used together.

The lowest SHAP importance is observed for γd. This does not imply that γd has no physical influence on MR. γd affects soil compactness, void ratio, and particle contact conditions. However, in the present dataset, γd exhibits a relatively narrow variation range compared with qc and w. In addition, part of the density-related stiffness information may already be reflected indirectly by qc and w. Consequently, the model attributes a smaller independent contribution to γd. Overall, the SHAP ranking is consistent with the geomechanical mechanisms controlling the MR of cohesive subgrade soils: in-situ mechanical resistance and moisture condition dominate the prediction, while fs and γd provide supplementary information.

4  Conclusions

This study developed an explainable GAN-Optuna-MLP framework for predicting the MR of soils under small-sample conditions. The proposed framework integrates GAN-based data augmentation, Optuna-based hyperparameter optimization, MLP-based nonlinear regression, and SHAP-based model interpretation. In addition, Optuna-RF and Optuna-XGBoost models were introduced as benchmark models to further evaluate the predictive performance of the proposed framework. Based on the results obtained in this study, the following conclusions can be drawn.

(1)   The GAN-generated samples showed good consistency with the original data in terms of distribution characteristics, descriptive statistics, and correlation patterns. This indicates that the GAN model can effectively capture the statistical characteristics of the limited original dataset and generate reasonable synthetic samples for model training.

(2)   Compared with the Optuna-MLP model without GAN-based data augmentation, the proposed GAN-Optuna-MLP model substantially improved the prediction accuracy. On the independent testing set, GAN-based data augmentation increased R2 by 6.9%, reduced MAE by 24.9%, and decreased MAPE by 25.1%. These improvements confirm that GAN-generated samples can effectively alleviate the data scarcity problem and enhance the generalization ability of the MLP model.

(3)   Compared with the optimized ensemble learning models, including Optuna-RF and Optuna-XGBoost, the proposed GAN-Optuna-MLP model achieved comparable predictive accuracy and slightly better error-based performance. This suggests that although RF and XGBoost are strong benchmark models for small tabular datasets, the proposed framework enables MLP to achieve competitive performance through data augmentation and hyperparameter optimization.

(4)   The SHAP analysis provided an interpretable understanding of the prediction mechanism of the developed model. The results showed that cone tip resistance qc was the most influential variable for MR prediction, while dry density γd had a relatively smaller contribution.

Overall, the proposed GAN-Optuna-MLP framework provides an effective and interpretable data-driven approach for MR prediction under limited-data conditions. Nevertheless, the trained model was developed based on a relatively small database of cohesive subgrade soils and is not directly claimed to be universally applicable to all soil types. Future studies should further validate and extend the proposed framework using larger multi-source databases, different soil categories, additional feature variables, and field-scale engineering cases.

Acknowledgement: Not applicable.

Funding Statement: This work was supported by the Postgraduate Research and Practice Innovation Program of Jiangsu Province (No. KYCX23-0686).

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Jian Liu and Yuedong Wu; methodology, Huiguo Wu; software, Huiguo Wu; validation, Jian Liu and Yuedong Wu; investigation, Huiguo Wu; resources, Yuedong Wu; data curation, Huiguo Wu; writing—original draft preparation, Huiguo Wu; writing—review and editing, Yuedong Wu; visualization, Huiguo Wu; supervision, Yuedong Wu; project administration: Jian Liu; funding acquisition, Yuedong Wu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The primary dataset used in this study was derived from the 124-sample clay database reported by Liu et al. [3], where the input variables and resilient modulus values used for modeling are explicitly listed. The code and processed modeling files generated during the current study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Salour F, Erlingsson S. Resilient modulus modelling of unsaturated subgrade soils: laboratory investigation of silty sand subgrade. Road Mater Pavement Des. 2015;16(3):553–68. doi:10.1080/14680629.2015.1021107. [Google Scholar] [CrossRef]

2. Zaman M, Solanki P, Ebrahimi A, White L. Neural network modeling of resilient modulus using routine subgrade soil properties. Int J Geomech. 2010;10(1):1–12. doi:10.1061/(asce)1532-3641(2010)10:1(1). [Google Scholar] [CrossRef]

3. Liu S, Zou H, Cai G, Bheemasetti TV, Puppala AJ, Lin J. Multivariate correlation among resilient modulus and cone penetration test parameters of cohesive subgrade soils. Eng Geol. 2016;209(2):128–42. doi:10.1016/j.enggeo.2016.05.018. [Google Scholar] [CrossRef]

4. Park HI, Kweon GC, Lee SR. Prediction of resilient modulus of granular subgrade soils and subbase materials using artificial neural network. Road Mater Pavement Des. 2009;10(3):647–65. doi:10.1080/14680629.2009.9690218. [Google Scholar] [CrossRef]

5. Salour F, Erlingsson S, Zapata CE. Modelling resilient modulus seasonal variation of silty sand subgrade soils with matric suction control. Can Geotech J. 2014;51(12):1413–22. doi:10.1139/cgj-2013-0484. [Google Scholar] [CrossRef]

6. Mehrotra A, Abu-Farsakh M, Gaspard K. Development of subgrade Mr constitutive models based on physical soil properties. Road Mater Pavement Des. 2018;19(1):56–70. doi:10.1080/14680629.2016.1235506. [Google Scholar] [CrossRef]

7. Sadrossadat E, Heidaripanah A, Ghorbani B. Towards application of linear genetic programming for indirect estimation of the resilient modulus of pavements subgrade soils. Road Mater Pavement Des. 2018;19(1):139–53. doi:10.1080/14680629.2016.1250665. [Google Scholar] [CrossRef]

8. Sadrossadat E, Ghorbani B, Zohourian B, Kaboutari M, Rahimzadeh Oskooei P. Predictive modelling of the MR of subgrade cohesive soils incorporating CPT-related parameters through a soft-computing approach. Road Mater Pavement Des. 2020;21(3):701–19. doi:10.1080/14680629.2018.1527241. [Google Scholar] [CrossRef]

9. Ikeagwuani CC, Nwonu DC, Nweke CC. Resilient modulus descriptive analysis and estimation for fine-grained soils using multivariate and machine learning methods. Int J Pavement Eng. 2022;23(10):3409–24. doi:10.1080/10298436.2021.1895993. [Google Scholar] [CrossRef]

10. Mousavi SH, Gabr MA, Borden RH. Resilient modulus prediction of soft low-plasticity Piedmont residual soil using dynamic cone penetrometer. J Rock Mech Geotech Eng. 2018;10(2):323–32. doi:10.1016/j.jrmge.2017.10.007. [Google Scholar] [CrossRef]

11. Pahno S, Yang JJ, Kim SS. Use of machine learning algorithms to predict subgrade resilient modulus. Infrastructures. 2021;6(6):78. doi:10.3390/infrastructures6060078. [Google Scholar] [CrossRef]

12. Kardani N, Aminpour M, Nouman Amjad Raja M, Kumar G, Bardhan A, Nazem M. Prediction of the resilient modulus of compacted subgrade soils using ensemble machine learning methods. Transp Geotech. 2022;36(1):100827. doi:10.1016/j.trgeo.2022.100827. [Google Scholar] [CrossRef]

13. Khawaja L, Asif U, Onyelowe K, Al Asmari AF, Khan D, Javed MF, et al. Development of machine learning models for forecasting the strength of resilient modulus of subgrade soil: genetic and artificial neural network approaches. Sci Rep. 2024;14(1):18244. doi:10.1038/s41598-024-69316-4. [Google Scholar] [CrossRef]

14. Bardhan A, Kumar Suman S, Kumar S, Lekhraj, Asteris PG. An efficient framework of optimized ensemble paradigm for estimating resilient modulus of subgrades. Transp Geotech. 2024;48:101315. doi:10.1016/j.trgeo.2024.101315. [Google Scholar] [CrossRef]

15. Zhang P, Yin ZY, Jin YF. Machine learning-based modelling of soil properties for geotechnical design: review, tool development and comparison. Arch Comput Meth Eng. 2022;29(2):1229–45. doi:10.1007/s11831-021-09615-5. [Google Scholar] [CrossRef]

16. Yaghoubi E, Yaghoubi E, Khamees A, Vakili AH. A systematic review and meta-analysis of artificial neural network, machine learning, deep learning, and ensemble learning approaches in field of geotechnical engineering. Neural Comput Appl. 2024;36(21):12655–99. doi:10.1007/s00521-024-09893-7. [Google Scholar] [CrossRef]

17. Hong S, Lee MH, Yoo BS, Kwak TY, Kim SR. Application of deep learning algorithms for predicting consolidation settlement. KSCE J Civ Eng. 2025;29(1):100072. doi:10.1016/j.kscej.2024.100072. [Google Scholar] [CrossRef]

18. Głuchowski A. Prediction of resilient modulus value of cohesive and non-cohesive soils using artificial neural network. Materials. 2024;17(21):5200. doi:10.3390/ma17215200. [Google Scholar] [CrossRef]

19. Raiaan MAK, Sakib S, Fahad NM, Al Mamun A, Rahman MA, Shatabda S, et al. A systematic review of hyperparameter optimization techniques in Convolutional Neural Networks. Decis Anal J. 2024;11(1):100470. doi:10.1016/j.dajour.2024.100470. [Google Scholar] [CrossRef]

20. Xu B, Yang G. Interpretability research of deep learning: a literature survey. Inf Fusion. 2025;115(11):102721. doi:10.1016/j.inffus.2024.102721. [Google Scholar] [CrossRef]

21. Jia F, Zhang J, Shen J, Wu L, Ma S. Compaction quality assessment of road subgrades using explainable deep graph learning framework. Comput Geotech. 2024;176(4):106795. doi:10.1016/j.compgeo.2024.106795. [Google Scholar] [CrossRef]

22. Goodfellow I, Pouget-Abadie J, Mirza M, Xu B, Warde-Farley D, Ozair S, et al. Generative adversarial networks. Commun ACM. 2020;63(11):139–44. doi:10.1145/3422622. [Google Scholar] [CrossRef]

23. Al-qerem A, Alsalman YS, Mansour K. Image generation using different models of generative adversarial network. In: Proceedings of the 2019 International Arab Conference on Information Technology (ACIT); 2019 Dec 3–5; Al Ain, United Arab Emirates. p. 241–5. doi:10.1109/ACIT47987.2019.8991120. [Google Scholar] [CrossRef]

24. Mukhiddin T, Lee W, Lee S, Rashid T. Research issues on generative adversarial networks and applications. In: Proceedings of the 2020 IEEE International Conference on Big Data and Smart Computing (BigComp); 2020 Feb 19–22; Busan, Republic of Korea. p. 487–8. doi:10.1109/bigcomp48618.2020.00-19. [Google Scholar] [CrossRef]

25. Jiang SH, Zhu GY, Wang ZZ, Huang ZT, Huang J. Data augmentation for CNN-based probabilistic slope stability analysis in spatially variable soils. Comput Geotech. 2023;160(2):105501. doi:10.1016/j.compgeo.2023.105501. [Google Scholar] [CrossRef]

26. Chen M, Kang X, Ma X. Deep learning—based enhancement of small sample liquefaction data. Int J Geomech. 2023;23(9):04023140. doi:10.1061/ijgnai.gmeng-8381. [Google Scholar] [CrossRef]

27. Akiba T, Sano S, Yanase T, Ohta T, Koyama M. Optuna: a next-generation hyperparameter optimization framework. In: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining; 2019 Aug 4–8; Anchoragek, AK, USA. p. 2623–31. doi:10.1145/3292500.3330701. [Google Scholar] [CrossRef]

28. Sheikh IR, Zhang M, Sun X, Chen C, Chen X, Dong Z, et al. Interpretable machine learning framework for resilient modulus estimation using LTPP data for pavements. Case Stud Constr Mater. 2025;23(5):e05403. doi:10.1016/j.cscm.2025.e05403. [Google Scholar] [CrossRef]

29. Chen Y, Kadkhodaei MH, Zhou J. Development of the Optuna-NGBoost-SHAP model for estimating ground settlement during tunnel excavation. Undergr Space. 2025;24(5):60–78. doi:10.1016/j.undsp.2025.03.006. [Google Scholar] [CrossRef]

30. Lundberg SM, Lee SI. A unified approach to interpreting model predictions. In: Guyon I, von Luxburg U, Bengio S, Wallach H, Fergus R, Vishwanathan S et al., editors. Advances in neural information processing systems 30. Long Beach, CA, USA: Curran Associates; 2017. p. 4765–74. [Google Scholar]

31. Li X, Liu W, Xu C, Liu N, Feng S, Zhang X, et al. Explainable machine learning-based prediction model for dynamic resilient modulus of subgrade soils. Transp Geotech. 2024;49(4):101415. doi:10.1016/j.trgeo.2024.101415. [Google Scholar] [CrossRef]