Advancements in predicting soil liquefaction susceptibility: a comprehensive analysis of ensemble and deep learning approaches

Semiempirical procedure

A semiempirical field-based procedure is used in this analysis to analyze the liquefied and non-liquefied case histories. The advantage of this approach is that it uses theoretical ideas and experimental results to establish the foundation for the analysis methodology and its constituent parts. The cyclic stress ratio (\({CSR}_{7.5}\)) was adjusted to a standard earthquake magnitude of 7.5, which helps in comparing the liquefaction potential across different seismic events. Equation (1) is used to determine \({CSR}_{7.5}\) at a depth z below the surface of the ground, as proposed by Seed and Idriss⁶. By including a magnitude scaling factor (\(MSF\)), the equivalent number of stress cycles is also adjusted for earthquakes of various magnitudes.

The factor 0.65 is a constant used in the liquefaction potential evaluation, representing an empirical reduction factor on the basis of observations and studies of earthquake-induced soil liquefaction. To standardize the cyclic stress ratio \((CSR)\) induced by an earthquake of magnitude \(M\) to an equivalent \(CSR\) for an earthquake with a magnitude of \(M=7.5\), a magnitude scaling factor \((MSF)\) is applied. The \(MSF\) factor adjusts the CSR for different earthquake magnitudes, as different magnitudes result in different durations and numbers of loading cycles. The MSF typically reduces the CSR for larger magnitudes. MSF is used to account for the duration effects in triggering soil liquefaction, specifically considering the number and relative amplitudes of loading cycles. This factor adjusts for the influence of earthquake magnitude on the potential for liquefaction by addressing the cumulative damage effects associated with longer-duration shaking. In this study, the MSF values were calculated via the methodology proposed by Idriss and Boulanger⁷⁶, which provides distinct MSF estimates for sandy and clayey soils. Their approach offers a nuanced assessment of liquefaction potential by incorporating soil type-specific characteristics and earthquake duration impacts, enhancing the accuracy of the liquefaction susceptibility evaluation, which is given by the following relationship presented in Eqs. (2) and (3):

However, the overburden correction factor was evaluated in terms of \({P}_{a}\) and \({\sigma {\prime}}_{vo}\) via the mathematical relationship presented in Eqs. (4) and (5).

It is necessary to use Eq. (4); however, this process can be simplified by selecting the automatic iteration option in an Excel spreadsheet. The number of q_c1ncs was limited to 21–254 for Eq. (5). The shear stress reduction factor can be estimated via the parameters of the earthquake magnitude (M) and depth (z). The stress reduction coefficient is evaluated via the following empirical formula presented in Eqs. (6) to (8):

However, Boulanger⁷⁷ established the basis for the \({K}_{\sigma }\) connection by demonstrating that the CRR for clean reconstituted sand in the laboratory may be related to the relative state parameter index of the sand. Idriss and Boulanger⁷⁶ suggested representing the ensuing \({K}_{\sigma }\) connection in terms of \({q}_{c1Ncs}\), as presented in Eqs. (9) and (10).

By limiting \({q}_{c1Ncs}\) to \(\le\) 211, the coefficient \({C}_{\sigma }\) can be reduced to its maximum value of 0.3. The equivalent clean sand adjustments for the fine content present in the soil were also considered in this study. According to current CPT-based research, the liquefaction case histories reveal that when the fines content (\(FC\)) increases, the liquefaction-triggering correlations move to the left. The equivalent clean sand adjustments, \(\Delta {q}_{c1N}\), are empirically derived from the liquefaction case history data. The adjusted expression for equivalent clean sand in the CPT is as follows in Eq. (11):

where \(FC\) denotes the percentage of fine content.

Finally, fine content and soil classification estimation analyses were performed by estimating the CPT tip resistance and sleeve friction ratio. The CPT tip resistance and sleeve friction ratio are functions of the soil behavior type index (\({I}_{C}\)), which is a function of \(FC\) and soil categorization. Robertson and Wride⁷⁸ suggested that the \({I}_{C}\) term be calculated via Eq. (12).

where Q and F are the normalized tip and sleeve friction ratios computed via Eqs. (13) and (14), respectively.

By first regressing \({I}_{C}\) against FC via the combined datasets to produce the least-squares fit, the connection for calculating \(FC\) was constructed via Eq. (15).

where \({C}_{FC}\) is a fitting parameter that can be modified depending on site-specific data.

The equation for the CPT-based data for calculating the cyclic resistance ratio (\(CRR\)) is as follows in Eq. (16).

where \({q}_{c1Ncs}\) is the equivalent clean sand adjustment given as follows in Eqs. (17) and (18).

where \({q}_{c1N}\) is the penetration resistance obtained from the same sand at an overburden stress of 1 atm when the other parameters remain constant.

Finally, the liquefaction safety factor is defined via Eq. (19).

Using Eq. (19), the factor of safety (FOS) against liquefaction was calculated for each dataset. A FOS less than 1 indicates that the soil is susceptible to liquefaction, suggesting a high risk of failure. Conversely, an FOS greater than 1 signifies that the soil is resistant to liquefaction, indicating stability under seismic conditions.

This study extends beyond conventional methodologies by exploring the potential of ensemble machine learning (ML) and deep learning (DL) techniques for predicting liquefaction. These data-driven approaches enable the extraction of intricate patterns from extensive datasets comprising soil properties and seismic parameters, potentially increasing the precision and automation of predictions. Figure 3 illustrates the methodology flowchart adopted in this investigation. A diverse array of algorithms is scrutinized, encompassing ensemble ML models such as XGBoost and RF models. Furthermore, this research delves into the capabilities of deep neural networks (DNNs), including LSTM and Bi-LSTM networks, for capturing prolonged dependencies. The efficacy of these algorithms will undergo meticulous assessment via a comprehensive dataset incorporating soil properties and seismic data.

eXtreme gradient boosting (XGBoost)

XGBoost XGBoost stands out as an advanced gradient boosting decision tree (GBDT) algorithm, which is an open-source library offering the application of ML algorithms for both regression and classification tasks in several fields^79,80. The library is recognized for its efficiency, flexibility, and portability, supporting multiple programming languages such as C +  + , Python, and R Studio. The XGBoost algorithm builds a sequence of classification or regression trees, known as CART, which serve as weak learners. These trees are sequentially combined to create the final prediction model. In line with other boosting techniques, XGBoost incrementally constructs regression trees, refines the model step-by-step and ensures that each tree minimizes the average loss function value across all steps in the training set. XGBoost enhances predictive accuracy by combining multiple weak learners, typically decision trees. Iteratively addressing previous prediction errors, each new model learns from residuals to improve overall accuracy. The inclusion of a regularization term in its objective function helps reduce overfitting and manage model complexity. The algorithm constructs a series of classification or regression trees (CARTs) step by step, and XGBoost minimizes the average loss function value at each step, forming a robust final predictive model.

Specifically, the model is trained on the dataset \(D=\left\{{x}_{i}, {y}_{i}\right), where i=\left(1, 2 ,3\dots n\right) and {x}_{i}\in {R}^{m}\text{ is the input vector having m number of input variables}\), and the output vector is denoted as \({y}_{i}\in R\). XGBoost uses the following objective function and regularization term presented in Eqs. (20) and (21):

where \(\Omega \left(.\right)\) and \(L(.)\) denote the regularization term and loss function, respectively. The loss function \(L(.)\) predicts the liquefiable and non-liquefiable conditions for a given training sample. In this context, “T” denotes the number of leaves in a decision tree, whereas \({\omega }_{j}\) represents the weight assigned to each leaf. By applying a second-order Taylor series expansion, Chen and Guestrin⁸¹ derived an approach to optimize this loss function, as illustrated in Eq. (22).

In this context, \(g_{i}\) and \(h_{i}\) represent the first and second derivatives of the loss function, respectively. The final predictive model is trained by incorporating y as the estimation for the \({i}^{th}\) instance at the \({t}^{th}\) iteration.

This equation assesses the suitability of a tree for the current step, but optimal values can be determined only after the tree structure is established. Owing to the impracticality of evaluating all possible structures, XGBoost constructs trees iteratively. XGBoost incorporates a regularization term \(\Omega\) and allows users to adjust two parameters: the maximum depth and the learning rate \(\eta\). The maximum depth, ranging from 0 to ∞ with a default of 8, limits tree depth, whereas (0< \(\eta\)  < 1) scales the prediction of each tree, reducing overfitting and enhancing model performance.

Random forest (RF)

The RF model is a popular machine learning algorithm utilized to solve several complex and nonlinear problems at various scales. In this study, the RF model was utilized to assess the liquefaction success of soil by analyzing various input features related to soil and seismic conditions. The RF is constructed with the help of several ensembles of decision trees (DTs), where each DT is trained via a random data sample and input features. The diversity in the training subsets helps capture complex, nonlinear relationships between input variables and the likelihood of liquefaction. In the context of soil liquefaction prediction, the random forest model can handle a variety of input parameters, such as soil and seismic parameters (e.g., variables D, FC, \({a}_{max}\), \({q}_{c}\), \({\sigma }_{v}\), \({\sigma {\prime}}_{v}\), \({M}_{w}\), and \({R}_{f}\)). By voting on the predictions from multiple decision trees, the random forest model provides a robust estimate of the liquefaction potential of the soil.

One of the key advantages of using the random forest model in this context is its ability to handle large datasets and manage missing data without significant performance loss. It also offers insights into the importance of different features in predicting liquefaction, which can be valuable for understanding the key factors influencing soil behavior during seismic events. A total of 336 data cases were used for training the RF model, and 143 cases were used for testing. Using the training dataset, the link between the input factors and liquefaction was modeled, and the accuracy of the predictions was then determined via the testing dataset. The RF model can handle large databases and runs efficiently on thousands of input variables⁸². Overall, the random forest model is a powerful tool for assessing liquefaction risk, assisting in the development of effective mitigation strategies.

Long short-term memory

The long short-term memory (LSTM) model is a specialized type of recurrent neural network (RNN) that excels in learning and predicting sequences of data LSTM networks and was introduced by Hochreiter and Schmidhuber⁸³. Since their introduction, LSTMs have been widely applied across diverse domains, including language modeling, speech recognition, and time series prediction. It is particularly useful in scenarios where temporal dependencies and sequential patterns are crucial, making it a promising approach for predicting the liquefaction potential of soil. In the context of soil liquefaction, the LSTM model can be used to analyze several soil and earthquake parameters that influence the liquefaction potential of soil during seismic events. Unlike traditional models, LSTMs can capture long-range dependencies and patterns in the data owing to their unique architecture, which includes memory cells and gating mechanisms. These components enable the model to retain relevant information from previous time steps and selectively forget irrelevant details, which is critical in understanding the cumulative effects of seismic events on soil stability. By leveraging historical data, such as past earthquake records and associated soil responses, the LSTM model can predict the likelihood of liquefaction under future seismic conditions. This predictive capability is particularly valuable for early warning systems and the planning of mitigation strategies. Moreover, LSTMs can incorporate various features, including soil properties, seismic characteristics, and environmental factors, providing a comprehensive assessment of liquefaction potential.

The memory cell, which is essential to long short-term memory (LSTM) networks, is engineered to maintain its state for extended periods of time. The network can then learn the dependencies that persist over time. Long short-term memory (LSTM) circuits regulate data entry and exit via three distinct types of gates: forget, input, and output. After deciding which input values to update the cell state with, it controls the output depending on the cell state and decides what information to discard. As relative information travels down the sequence chain, the cell state serves as a transport highway. It follows the entire chain in a straight line, with only a few small linear exchanges, enabling data to be transferred across numerous time steps.

The forget gate in an LSTM network is responsible for determining which information should be retained and which should be discarded. It processes the previous hidden state, \({h}_{t-1}\), and the current input, \({x}_{t}\), through a sigmoid activation function, as described in Eq. (23). The resulting output is a vector of values between 0 and 1, corresponding to each element in the cell state \({c}_{t-1}\)​. A value closer to 1 indicates that the information should be retained, whereas a value closer to 0 suggests that the information can be forgotten, thus enabling the network to manage long-term dependencies effectively.

where \(f\left(t\right)\) represents the output of the forget gate and where \({P}_{f}\) and \({b}_{f}\) represent the weight matrix and bias vector of the forget gates from the input layer, respectively.

The input gate decides which values will be updated. It has two parts: a sigmoid layer (deciding which values to update) and a \(tanh\) layer (creating a vector of new candidate values) via Eqs. (24) and (25).

The old cell state \({C}_{t-1}\) is updated into the new cell state \({C}_{t}\). This involves forgetting some parts of the old cell state and adding new candidate values via Eq. (26).

The output gate determines what the next hidden state \({h}_{t}\) should be. This is based on the cell state and involves passing it through a tanh function and then multiplying it by the output gate.

LSTM networks are powerful tools for handling sequences of data and have proven effective in a wide range of applications because of their ability to maintain long-term dependencies and manage the vanishing gradient problem typical of standard RNNs.

Bidirectional long short-term memory (Bi-LSTM)

Bidirectional long short-term memory (BI-LSTM) networks represent advanced iterations of traditional LSTM networks engineered to enhance the comprehension of input sequences. Unlike standard LSTM networks, which process data in a unidirectional manner, BI-LSTM networks analyze sequences in both forward and backward directions. This bidirectional processing enables the model to access contextual information from both past and future inputs simultaneously. This architecture is especially advantageous in tasks that necessitate a comprehensive understanding of the full context, as it allows the network to capture dependencies and nuances that may otherwise be overlooked⁸³.

In a Bi-LSTM network, each input sequence undergoes dual processing through two distinct LSTM networks. The first, known as the forward LSTM, processes the sequence from the initial element to the final element (from left to right). Conversely, the second network, the backward LSTM, processes the sequence in the reverse order, from the final element to the initial one (right to left). At each time step, the outputs generated by the forward and backward LSTM networks are concatenated. This concatenation allows the BI-LSTM to leverage information from both preceding and succeeding elements, thereby enhancing the model’s contextual understanding. This concatenation provides a comprehensive representation of the input at each time step, combining information from both directions. The input sequence is provided by \(X=[{x}_{1},{x}_{2},…,{x}_{T}]\), and the forward LSTM processes it in the original order to produce the forward hidden states \({\overrightarrow{h}}_{t}\)​ presented in Eq. (29).

After the forward LSTM process is forwarded, the reverse sequence is reserved, and the backward LSTM process produces the backward hidden states \({\overleftarrow{h}}_{t}\) via Eq. (30).

At each time step \(t\), the forward and backward hidden states are concatenated to form the final hidden state \({h}_{t}\), as presented in Eq. (31).

The concatenated hidden states \({h}_{t}\) are then used for subsequent layers or the final output layer, depending on the specific task. In this work, \({h}_{t}\) are used for the liquefaction classification system.

By processing the sequence in both directions, BI-LSTM can capture information from the entire sequence, leading to a better understanding of the context to overcome the limitations of the LSTM model. BI-LSTMs often achieve better performance on tasks that require context from both past and future inputs. BI-LSTM networks improve upon traditional LSTMs by integrating information from both preceding and succeeding contexts. This dual-directional approach makes them exceptionally effective for tasks requiring comprehensive sequence understanding, such as language processing and time series analysis. Consequently, BI-LSTMs provide richer, more accurate representations of input data, enhancing overall performance and contextual accuracy.

Performance assessment

The performance assessment of classification machine learning models is a critical aspect of the model development process. It involves evaluating how well a model predicts the class labels for new, unseen data. This assessment helps in understanding the model’s accuracy, robustness, and generalizability. A confusion matrix is a table that helps visualize the performance of a classification algorithm. It compares the actual target values with those predicted by the model. The calculation of several performance metrics is performed via the two-class confusion matrix shown in Fig. 4. True positives (TPs) are examples of the number of correctly predicted liquefiable cases that the algorithm properly recognizes as positive, whereas the negative components that are accurately classified as negative are known as true negatives (TNs). On the other hand, false negatives (FNs) are cases of positive data for which the ML models incorrectly classify the data as negative. The positive components that are incorrectly forecasted as negatives are known as false positives (FPs).

Various performance fitness error matrices (PFEMs) were calculated to assess the accuracy and reliability of the proposed models. The PFEMs included in this study are accuracy, recall, specificity, precision, F1 score, MCC, and balance accuracy (BA). The expressions for these PFEMs are presented in Eqs. (32) to (38).

ML model hyperparameter configuration

GridSearchCV is a powerful tool in the sci-kit-learn library that is used for hyperparameter tuning. This allows us to exhaustively search through a specified parameter grid and determine the optimal parameters for the XGBoost and RF algorithms. When the XGBoost and RF algorithms are used, GridSearchCV helps in finding the best combination of hyperparameters that leads to the highest performance of the model. Using GridSearchCV effectively requires a good understanding of the model and the problem at hand, allowing us to select an appropriate range of hyperparameters to explore. By systematically testing combinations of parameters such as the number of estimators, maximum tree depth, and learning rate, GridSearchCV identifies the configuration that yields the best model performance on the validation data. This approach ensures a more objective and data-driven selection of hyperparameters, reducing the reliance on manual trial-and-error. As shown in Tables 2 and 3, the optimal values were determined based on defined search ranges for each algorithm. For example, the XGBoost model achieved optimal performance with 150 trees, a maximum depth of 6, and a learning rate of 0.1, whereas the RF model performed best with 200 trees and a maximum depth of 6. Although GridSearchCV is primarily applicable to traditional ML models, the hyperparameters for deep learning models such as LSTM and BI-LSTM were selected through empirical experimentation, as detailed in Table 4. These included settings such as the number of LSTM units, dropout rate, batch size, and training epochs. Together, these tuning strategies significantly contribute to enhancing model accuracy and robustness in the assessment of liquefaction potential.