An Application of Particle Filter for Parameter Estimation and Prediction in Geotechnical Engineering

Abstract

This chapter presents an application example of a nonlinear Kalman Filters (KFs), i.e., Particle Filter (PF), for state (or parameters) estimation and prediction of a dynamical system in geotechnical engineering. First key characteristics of dynamical systems in geotechnics, which need to be considered in filtering, are described by showing some figures, and why PF is necessary for geotechnical applications is explained. Then, a detailed algorithm and implementation of PF for geotechnical problems are presented with key equations. The PF is demonstrated through a case history focusing on deformation behavior of a ground due to embankment construction. The PF is applied to estimation of geotechnical parameters and predictions of future settlement behavior of the ground to discuss the effectiveness of the PF in geotechnical engineering. The results of the case history have shown that PF has presented great promise as an accurate parameter identification for a nonlinear dynamic model. The simulation with the identified parameters predicts the actual measurement data with high accuracy even though a limited amount of measurement data was used in identification stage. The PF provides more information on estimates than optimization methods because the estimates are obtained in the form of probability density functions (PDFs). This characteristic can contribute to risk analysis and reliability-based decision-making in geotechnical practice.

1. Introduction

Filtering refers to the process of estimating state or parameters and prediction of dynamic systems based on noisy monitoring data, and it has significant applications in many different fields. There is no exception in geotechnical engineering: there are many applications of filtering for geotechnical problems, particularly for enhancing the accuracy and reliability of numerical simulations. Filtering helps in fine-tuning numerical simulation models by continuously updating them as new data becomes available. This iterative process enhances the accuracy and reliability of the models, leading to better decision-making in engineering projects. One of the key benefits of filtering is its ability to quantify uncertainties in model predictions. Geotechnical engineering involves dealing with uncertain soil properties, geological variations, and environmental conditions [1]. Filtering provides a probabilistic framework to express these uncertainties, which is crucial for risk assessment and management in engineering projects.

One of the most well-known methods or algorithms to achieve filtering is KF. In 1960, R.E. Kalman published a famous paper describing a recursive solution to the discrete-data linear filtering problem [2]. This algorithm has been applied to the Apollo program. Since that time, due in large part to advances in digital computing, the KF has been the subject of extensive research and application, particularly in the area of time-series analysis. The KF is a set of mathematical equations that provides an efficient computational recursive solution of the least-squares method (LSM). The filter is very powerful in several aspects: it supports estimations of past, present, and even future states, and it can do so even when the precise nature of the modeled system is unknown.

KF also received significant attention in geotechnical engineering around 1980s. Murakami and Hasegawa [3] applied KF to a geotechnical one-dimensional time-series problem (consolidation settlement problem) for updating geotechnical parameters and predicting consolidation settlement of soft clay ground. Murakami and Hasegawa [4, 5] extended their method to two-dimensional settings: they integrated KF and Finite Element Method (FEM), which is known as “Kalman Filter FEM”, for solving broad practical/realistic time-series estimation problems in geotechnical engineering. In addition, they utilized it successfully to solve the optimal location of observation points [6, 7]. After these pioneering works, several researches combining KF with numerical simulation models were reported in geotechnical engineering [8, 9, 10, 11].

Although the KF is a powerful tool for state (or parameter) estimations and prediction, it has some limitations which are as follows:

1. Linearity: Standard KF assumes that the system dynamics and measurement models are linear. This can be a significant limitation for systems that are inherently nonlinear.

2. Gaussianity: KF assumes that the process and measurement noise are both Gaussian. If the noise characteristics deviate significantly from a Gaussian distribution, the filter’s performance may degrade.

3. Robustness: KF may not be robust to outliers or sudden changes in the system dynamics, which can lead to divergence or inaccurate estimates.

Unfortunately, the characteristics of geotechnical problems can make the application of the KF challenging. First, many geotechnical systems exhibit nonlinear behavior, such as soil-structure interaction and material nonlinearity (e.g., plasticity, strain-hardening). Second, geotechnical systems often encounter non-Gaussian noise due to complex subsurface conditions, measurement errors, and model uncertainties. KF assumes Gaussian noise, which may not be suitable for accurately capturing the uncertainties in geotechnical problems.

Despite these challenges, KF and its variants are still valuable tools in geotechnical engineering for applications such as deformation prediction, slope stability analysis, and foundation performance assessment. The variants include Extended Kalman Filter (EKF) [12] for nonlinear systems, the Unscented Kalman Filter (UKF) [13], Ensemble Kalman Filter (EnKF) [14] and Particle Filter (PF) [15, 16, 17] for better handling of nonlinearities, and adaptive versions for systems with changing dynamics. Due to the increased need for nonlinear filters in geotechnical engineering, several attempts have been made to apply these nonlinear filters for enhancing the accuracy and reliability of numerical simulations, managing uncertainties, adapting to changing conditions, and supporting cost-effective and safe engineering solutions. Its integration into geotechnical practices represents a significant advancement in the field, enabling more informed and confident decision-making in complex engineering projects [18, 19, 20, 21, 22, 23].

This chapter describes some key characteristics of geotechnical problems that need to be considered in filtering and presents an application example of a variant of nonlinear KF for parameter estimation and predictions in geotechnical engineering. To the best of the authors’ knowledge, PF is the most suitable method for geotechnical applications because of its versatility for nonlinear and non-Gaussian problems. We believe that this chapter is beneficial not only for geotechnical engineers and researchers but also for those who learn standard Kalman filtering and other nonlinear Kalman filtering. This chapter is organized as follows: Section 2 describes the importance and the roles of state estimations and predictions based on monitoring data in geotechnical engineering considering characteristics of geotechnical problems. Section 4 shows an application example of PF in a numerical simulation of a typical geotechnical problem: predicting deformation behavior of a ground due to embankment construction. Section 5 briefly summarizes this chapter.

2. Characteristics of geotechnical problems

Understanding heterogeneity of grounds and variability of soil parameters is crucial in geotechnical engineering. This section outlines the importance of state estimation and prediction based on monitoring data by considering heterogeneity of grounds and uncertainties in geotechnical parameters.

The concept of ground heterogeneity is a fundamental aspect of geotechnical engineering. It refers to the variation in the composition, structure, and properties of the ground, including soil and rock layers. This heterogeneity can occur on various scales and has significant implications for construction, environmental projects, and geological studies. At least two types of heterogeneity in the context of geotechnical engineering, which are a) soil composition: variation in the types of soil (clay, sand, silt, and gravel) and their mixtures, and b) rock structure: differences in rock types, stratification, faulting, folding, and weathering. Figure 1 shows an example of heterogeneity of the ground created by a coupled Markov chain (CMC) method [24]. The heterogeneity can be due to natural processes like sedimentation, erosion, and weathering, or human activities such as excavation or filling. Therefore, the heterogeneity of the ground can be observed not only in natural grounds, but also artificial grounds such as earthen embankments for river dikes.

Figure 1.

An example of heterogeneity of a ground [24].

Variability (or uncertainty) of soil parameters is another fundamental aspect in geotechnical engineering that influences the design and analysis of geotechnical structures. These uncertainties arise due to the complex and variable nature of soil as a material, and they have significant implications for the safety, reliability, and cost-effectiveness of geotechnical projects. There are many potential sources of natural variability of soils, which are a) natural variability: soils are formed through long geological processes, leading to inherent variability in properties like grain size, mineral composition, density, strength, and permeability, b) sampling and testing limitations: soil properties are typically determined through field sampling and laboratory testing, but samples may not represent the entire site, and tests can have inherent errors or limitations, C) spatial variability: soil properties can vary significantly over short distances, both horizontally and vertically, making it challenging to obtain a comprehensive understanding of the subsurface conditions, and d) temporal changes: soil properties can change over time due to factors like weathering, organic matter decomposition, moisture content changes, and human activities. Figure 2 shows examples of variability of rock parameters. We collected compressive strength and elastic modulus of granite from several technical reports published in Japan. It is clear that the parameters show large variations despite all the data being from Japanese granite.

Figure 2.

Variability of rock properties.

The uncertainties dealt with in geotechnical engineering are classified as being either aleatory – the inherent variation given sufficient samples of the stochastic process, which can be characterized via a probability distribution, or epistemic – where there is insufficient information concerning the parameters of interest. The aleatory uncertainty is also called stochastic uncertainty or randomness. On the other hand, the epistemic uncertainty is referred to as knowledge uncertainty or reducible uncertainty. The categories of uncertainty in soil properties are shown in Figure 3. These uncertainties need to be considered in geotechnical engineering.

Figure 3.

Categories of uncertainty in soil properties [25].

The heterogeneity and the parameter variability make numerical simulations in geotechnical engineering challenging. As shown in Figure 1, natural (and even artificial) grounds are essentially heterogeneous, and it is challenging to make numerical models considering the actual heterogeneity in simulations. Most of practical simulations assume that the ground is homogeneous for simplicity. In addition, geotechnical researchers and engineers usually assume that the input parameters are “deterministic” despite actual parameters being quite uncertain (Figure 2). Geotechnical engineers select one parameter from the dataset of the parameters and tend to choose “conservative” parameters to avoid the worst-case scenario.

Rather than using single, deterministic values for soil parameters and ground conditions, probabilistic approaches are necessary to deal with the uncertainties. In general, uncertainty quantification is necessary, as current real-world data is insufficient and incomplete. A renowned statistician C. R. Rao stated that the logical equation is important in decision-making under uncertainties (Figure 4). If we have to make decisions under uncertainty, we cannot avoid mistakes. If mistakes cannot be avoided, we had better know how often we make mistakes [26]. Therefore, it is natural to consider the probabilistic methods as suitable for geotechnical applications.

Figure 4.

A logical equation stated by C.R. Rao [26].

As mentioned in the previous section, one of the key benefits of Bayesian updating is its ability to quantify uncertainties. Uncertainties in geotechnical engineering can be quantitatively modeled with PDFs or statistical models, and the initial probability or statistical model can be updated (improved) based on observation data. Bayesian updating provides a reasonable framework for this process.

Mechanical behavior of soils shows strong nonlinearity. The nonlinearity refers to the characteristic that the response of soil to stress or strain is not directly proportional. In other words, the relationship between the applied stress and the resulting strain in soil is not a straight line, but rather it exhibits a more complex, often curved relationship. Figure 5 shows an example of stress–strain relationship of a clayey soil. At the beginning of the curve, the relationship between stress and strain is linear. In this region, the soil deforms proportionally to the applied stress, and the deformation is usually reversible if the stress is removed. As the stress increases, the curve reaches a point where the soil begins to yield, meaning it starts to undergo plastic deformation. Beyond this point, the relationship between stress and strain is no longer linear. After the yield point, the curve enters the plastic region, where the soil continues to deform under constant or slightly increasing stress. If the stress continues to increase, the curve eventually reaches a peak. The failure point is characterized by a significant increase in strain with little or no increase in stress. Beyond the failure point, the stress decreases while the strain continues to increase.

Figure 5.

Example of nonlinear behavior of soil (stress–strain curve).

It is clear that the stress–strain curve in Figure 5 shows strong nonlinearlity which is a critical factor in geotechnical engineering that can affect the design and analysis of geotechnical structures such as foundations, slopes, retaining walls, and other soil-structure interaction problems. It requires sophisticated nonlinear models and analysis techniques to accurately predict soil behavior under various loading conditions.

Table 1 summarizes several filtering algorithms. As mentioned, standard KF is designed for linear and Gaussian problems and is difficult to apply for practical geotechnical problems. Although nonlinear filtering such as EKF, UKF, and EnKF can analyze nonlinear problems, they cannot capture complex non-Gaussian probability functions. As shown in Figure 2, histograms of soil and rock parameters often show non-Gaussian distributions. Hence, PF is the most preferable algorithm of filtering for geotechnical applications.

Algorithm	Problem or model	Likelihood (Observation noise)	Prior probability
Kalman filter (KF)	Linear	Gaussian	Gaussian
Extended Kalman filter (EKF)	Nonlinear	Gaussian	Gaussian
Unscented Kalman filter (UKF)	Nonlinear	Gaussian	Gaussian
Ensemble Kalman filter (EnKF)	Nonlinear	Gaussian	Gaussian
Particle filter (PF)	Nonlinear	Non-Gaussian	Non-Gaussian

Table 1.

Some algorithms for sequential Bayesian filtering [7].

3. A versatile nonlinear filter: particle filter

As discussed in Section 2, geotechnical problems are fundamentally nonlinear and non-Gaussian, hence it is natural to use PF for geotechnical problems (see Table 1). This section provides a thorough description of the PF.

We consider a nonlinear and non-Gaussian state space model, which is represented in the following [16]:

$$x_t=f_t\left(x_{t-1}\right)+v_t$$E1

$$y_t=h_t\left(x_t\right)+w_t$$E2

where x_t = (x₁, x₂, …, x_k) and y_t = (y₁, y₂, …, y_l) indicate the state (or parameter) vector and the observation (or data) vector, and the subscripts k and l indicate the size of the vectors. f_t and h_t are the nonlinear operator (or function) and the observation operator. Subscript t indicates the discrete time. The vectors v_t and w_t are the vectors for system and observation errors respectively. The errors follow the normal (or Gaussian) distributions with zero mean which are defined by:

$$v_t\sim N\left(0,Q_t\right)$$E3

$$w_t\sim N\left(0,R_t\right)$$E4

where Q_t and R_t indicate covariance matrices.

This chapter focuses on the mechanical behavior of geomaterials which are commonly modeled by soil–water coupled FEM. When the soil–water coupled FEM is used, Eq. (1) is given by the following equation:

$$\left\{\begin{array}{c}{u}_t\\ p{}_t\end{array}\right\}=\left\{\begin{array}{c}{u}_{t-1}\\ p{}_{t-1}\end{array}\right\}+{\left[\begin{array}{cc}\left[\mathrm{K}\right]& {\left[{\mathrm{K}}_{\mathrm{v}}\right]}^{\mathrm{T}}\\ \left[{\mathrm{K}}_{\mathrm{v}}\right]& -\mathrm{\theta }\mathrm{\Delta }t\left[{\mathrm{K}}_{\mathrm{h}}\right]\end{array}\right]}^{-1}\left\{\begin{array}{c}\left\{\mathrm{F}\right\}+{\left[{\mathrm{K}}_{\mathrm{v}}\right]}^{\mathrm{T}}{p}_t\\ \left\{\mathrm{Q}\right\}+\left(1-\mathrm{\theta }\right)\left[{\mathrm{K}}_{\mathrm{h}}\right]p_t\end{array}\right\}+\left\{\begin{array}{c}{v}_t^{\mathrm{u}}\\ v_t^{\mathrm{p}}\end{array}\right\}$$E5

where

[K]: the tangent stiffness matrix defined by constitutive models

[K_v]: the matrix transforming the increment in nodal displacements to the volume change of each element

[K_h]: the matrix for seepage flow in soils

{F}: a vector representing applied force

{u_t}: the displacement vector

{p_t}: the ore pressure vector

θ: coefficient (0 < θ < 1)

{Q}: the vector representing increment in the volume rate of the water flow

{v_t^u}: the system error vector for {u_t}

{v_t^p}: the system error vector for {p_t}

The state (or parameter) vector x_t is defined by

$$x_t=\left\{\begin{array}{c}{u}_t\\ p_t\end{array}\right\}$$E6

If state variables u₁, u₂, and u₃ in {u_t} are directly observed, nonlinear function h_t can be written in matrix form as

$$H_t=\left[\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ 0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\end{array}\right]$$E7

We can describe y_t = H_tx_t using Eq. (7).

The recursive formula of probability densities for prediction and filtering is defined as follows:

Prediction (time update):

$$\begin{array}{lll}p\left(x_t\left|y_{1:t-1}\right.\right)& =& {\int }_{-\infty }^{\infty }p\left(x_t,x_t\left|y_{1:t-1}\right.\right){\mathit{dx}}_{t-1}={\int }_{-\infty }^{\infty }p\left(x_t\left|x_{t-1},y_{1:t-1}\right.\right)p\left(x_{t-1}\left|y_{1:t-1}\right.\right){\mathit{dx}}_{t-1}\\ & =& {\int }_{-\infty }^{\infty }p\left(x_t\left|x_{t-1}\right.\right)p\left(x_{t-1}\left|y_{1:t-1}\right.\right){\mathit{dx}}_{t-1}\end{array}$$E8

Filtering (observation update):

$$p\left(x_t\left|y_{1:t}\right.\right)=\frac{p\left(x_t,y_t\left|y_{1:t-1}\right.\right)}{p\left(y_t\left|y_{1:t-1}\right.\right)}=\frac{p\left(y_t\left|x_t\right.,y_{1:t-1}\right)\cdot p\left(x_t\left|y_{1:t-1}\right.\right)}{p\left(y_t\left|y_{1:t-1}\right.\right)}=\frac{p\left(y_t\left|x_t\right.\right)\cdot p\left(x_t,\left|y_{1:t-1}\right.\right)}{p\left(y_t\left|y_{1:t-1}\right.\right)}$$E9

where p(x_t|x_t-1) is the density of x_t given previous state vector x_t-1, p(y_t|x_t) is the density of y_t given x_t.

In linear and Gaussian problems, the KF provides updating recursions using Eqs (8) and (9). However, integral computations in (8) and (9) are usually intractable for practical problems which usually show nonlinearity and non-Gaussianity. This is the main reason to use numerical methods such as PF in geotechnical problems.

The PF expresses PDFs with a set of Monte Carlo (MC) samples called an “ensemble” that has weights, and each MC sample is called a “particle.” A filtered distribution at time t – 1, p(x_t-1|y_1:t-1), where y_1:t-1 denotes {y₁, y₂, …, y_t-1}, is redefined with “particles” and weights:

$$p\left(x_{t-1}\left|y_{1:t-1}\right.\right)\approx \frac{1}{N}\sum _{i=1}^N{w}_{t-1}^{\left(i\right)}\delta \left(x_{t-1}-x_{t-1\left|t-1\right.}^{\left(i\right)}\right)$$E10

where N is the number of particles and δ is the Dirac delta function. Using this particle or MC approximation, integral computation in (8) and (9) become summations.

The particle approximation for the prediction distribution p(x_t|y_1:t-1) at time t is given by the following equation:

$$\begin{array}{lll}p\left(x_t\left|y_{1:t-1}\right.\right)& =& {\int }_{-\infty }^{\infty }p\left(x_t\left|x_{t-1}\right.\right)p\left(x_{t-1}\left|y_{1:t-1}\right.\right){\mathit{dx}}_{t-1}\approx \sum _{i=1}^N{\int }_{-\infty }^{\infty }{w}_{t-1}^{\left(i\right)}\delta \left(x_{t-1}-x_{t-1\left|t-1\right.}^{\left(i\right)}\right)p\left(x_t\left|x_{t-1}\right.\right){\mathit{dx}}_{t-1}\\ & =& \sum _{i=1}^N{w}_{t-1}^{\left(i\right)}\delta \left(x_t-f_t\left(x_{t-1\left|t-1\right.}^{\left(i\right)},v_t^{\left(i\right)}\right)\right)=\sum _{i=1}^N{w}_{t-1}^{\left(i\right)}\delta \left(x_t-x_{t\left|t-1\right.}^{\left(i\right)}\right)\end{array}$$E11

where, {v_t⁽ⁱ⁾}_i=1^N is an independent and identically distributed (i.i.d.) sample set for Eq. (3). The calculation means that each particle for the prediction ensemble, x_t|t-1⁽ⁱ⁾, is obtained by the direct calculation of f_t(x_t-1|t-1⁽ⁱ⁾| v_t⁽ⁱ⁾).

Filtering

We obtain the ensemble approximation for filtered distribution, from Eqs (9) and (11), and observation y_t by the following calculation:

$$\begin{array}{lll}p\left(x_t\left|y_{1:t}\right.\right)& =& \frac{p\left(y_t\left|x_t\right.\right)\cdot p\left(x_t\left|y_{1:t-1}\right.\right)}{p\left(y_t\left|y_{1:t-1}\right.\right)}=\frac{p\left(y_t\left|x_t\right.\right)\cdot p\left(x_t\left|y_{1:t-1}\right.\right)}{{\int }_{-\infty }^{\infty }p\left(y_t\left|x_t\right.\right)\cdot p\left(x_t\left|y_{1:t-1}\right.\right){\mathit{dx}}_t}\\ & =& \frac{1}{{\sum }_{j}p\left(\left.y_t\right|x_{t\left|t-1\right.}^{\left(j\right)}\right)w_{t-1}^{\left(j\right)}}\sum _{i=1}^{N}p\left(\left.y_t\right|x_{t\left|t-1\right.}^{\left(i\right)}\right)w_{t-1}^{\left(i\right)}\delta \left(x_t-x_{t\left|t-1\right.}^{\left(i\right)}\right)\\ & =& \sum _{i=1}^N{\tilde{w}}_t^{\left(i\right)}{w}_{t-1}^{\left(i\right)}\delta \left(x_t-x_{t\left|t-1\right.}^{\left(i\right)}\right)=\sum _{i=1}^N{w}_t^{\left(i\right)}\delta \left(x_t-x_{t\left|t-1\right.}^{\left(i\right)}\right)\end{array}$$E12

where

$${\tilde{w}}_t^{\left(i\right)}=\frac{p\left(\left.y_t\right|x_{t\left|t-1\right.}^{\left(i\right)}\right)}{{\sum }_{j}p\left(\left.y_t\right|x_{t\left|t-1\right.}^{\left(j\right)}\right)w_{t-1}^{\left(j\right)}}$$E13

If the observation system is linear, p(y_t|x_t|t-1⁽ⁱ⁾) is given by Eq. (14):

$$p\left(\left.y_t\right|x_{\left.t\right|t-1}^{\left(i\right)}\right)=\frac{1}{{\left(2\pi \right)}^{m/2}\left|R_t\right|}exp\left[-\frac{{\left(y_t-H_t\left(x_{\left.t\right|t-1}^{\left(i\right)}\right)\right)}^{\mathrm{T}}{R}_t^{-1}\left(y_t-H_t\left(x_{\left.t\right|t-1}^{\left(i\right)}\right)\right)}{2}\right]$$E14

The rest of the problem is how to sample from p(x_t|y_1:t) which is called “sampling.” There are two ways of sampling methods which are Sequential Importance Resampling (SIR) and Sequential Importance Sampling (SIS). While classic PF algorithms use SIR [15, 16], SIS is known as a generalization of PF [17]. The pseudo-code of the two sampling methods is as follows:

Sequential Importance Resampling (SIR)

Step 1. Initial setup:

Generate an ensemble (set of particles) {x₀⁽¹⁾, x₀⁽²⁾,…, x₀^(N)} from the initial prior probability distribution p₀. Set t = 1.

Step 2. Prediction (Forecasting):

Each ensemble member or particle x_t-1⁽ⁱ⁾ evolves via nonlinear operator f_t given by Eq. (1).

Step 3. Filtering:

After measurement data y_t, is obtained, weight w_t⁽ⁱ⁾ is calculated.

Step 4. Resampling (based on filtered distribution):

Resample new particles {x_t⁽¹⁾, x_t⁽²⁾,…, x_t^(N)} N times from the set of particles or parameters x_t-1⁽ⁱ⁾. The set of samples approximates filtered distribution p(x_t|y_1:t).

Set t = t + 1 and go back to Step 2.

Sequential Importance Sampling (SIS)

Step 1. Initial setup:

Generate a set of particles {x₀⁽¹⁾, x₀⁽²⁾,…, x₀^(N)} from the initial (or prior) probability distribution.

Step 2. Prediction (Forecasting):

Each sample x_t-1⁽ⁱ⁾ evolves via nonlinear operator f_t given by Eq. (1)

Step 3. Filtering:

After measurement data y_t, is obtained, weight w_t⁽ⁱ⁾ is calculated.

Step 4. Weight (or probability) updating:

The set of weighted particles {x_t⁽ⁱ⁾} approximates filtered distribution p(x_t|y_1:t).

Set t = t + 1 and go back to Step 2.

Figure 6 compares the PDFs approximated using SIR and SIS.

Figure 6.

Filtered PDF approximated by particles.

The main question here is which sampling method, SIR or SIS, should be used for geotechnical applications. To answer the question, let us focus on the mechanical soil behavior. Soil undergoes both elastic and plastic deformation when load is applied to and reaches “critical state” under large shear deformation at constant volume and constant shear/normal stress conditions [27]. This mechanical behavior of soil can be well explained by the critical state constitutive models. In critical state soil mechanics, the mechanical state of soil is described by the stress parameters q and p and the specific volume (1 + e). The history of consolidation, which causes long-term settlement of the ground, is described by overconsolidation ratio p_m’/p’ (p_m’ is the stress at the intersection of the current swelling line with the normal consolidation line, Figure 7). In many theories for stress–strain behavior, the value of p_m’ defines the size of the state boundary surface in Figure 7 and it separates the elastic states inside from the elastoplastic states.

Figure 7.

Current stress and consolidation history (Modified from [28]).

The mechanical behavior or stress–strain response of soil not only depends on the current stress state, but also on the stress history that soil has undergone [28]. The schematic on the effect of the stress history is shown in Figure 8. Soil samples brought to the same initial states of q_i and p_i at zero (O), along different paths of CO and DO, are then loaded along the same path, OA. Figure 8b illustrates the stress–strain curves for the same loading path OA. Soil offers different stress–strain curves depending on the stress histories undergoes before OA loading, and this is commonly observed in laboratory tests for soil samples.

Figure 8.

Effect of recent stress history on current stiffness. (a) Stress paths of a soil, (b) stress-strain curves for the same loading path OA. (Modified from [28]).

Thus the stress or loading history is an essential factor for evaluating the mechanical (or stress–strain) behavior of soils. Let us refocus on the sampling methods herein. While the SIR algorithm abandons the information of stress history the soils have undergone due to resampling, the SIS keeps storing the information of stress history during the entire process of filtering. This implies that the SIS has natural ability to evaluate the mechanical behavior of soils considering stress history and it has high potential of Bayesian updating in geotechnical problems.

Selection of parameters to be identified and observed quantities [7]: The selection of parameters to be identified varies with the actual types of problems. For example, when a soft clay ground is improved by sand compaction piles or vacuum consolidation, the mechanical properties of the clay layer change after the ground improvement. Hence, the new, but unknown, properties of the improved layer, such as the hydraulic conductivity and compression index, can be a target for identification. Material constants that are difficult to measure, or parameters that significantly affect the future behavior of structures, are usually worthy of identification. Observed quantities, such as the settlement, lateral displacement, and pore water pressure, for such inverse analyses, must be sufficiently sensitive to the parameters to be identified. In practice, the measurable physical quantities are limited. Hence, we need to obtain and select the observed quantities that are relevant to the parameters to be identified.

Effective sample size [7]: When PF is used for Bayesian updating, attention must be paid to “degeneracy.” The term “degeneracy” refers to the scenario whereby, after a few iterations of the particle filter algorithm, only a small number of particles have significant weight. This means that only a few particles account for most of the posterior probability mass, making the approximation less accurate. Degeneracy happens when the weights of the particles become too unequal. There is a useful metric, which is called the effective sample size [29, 30], for discussing how severe the degeneracy is. The N_eff is defined by the following equation:

$$N_{\mathrm{eff}}=\frac{1}{{\sum }_{i=1}^N{\left(w_k^{\left(i\right)}\right)}^2}$$E15

where N is the total number of particles and w_k⁽ⁱ⁾ is the weight of the ith particle at the kth filtering step. The N_eff takes values from 1 (only one particle has significant weight) to N (all the particles have the same weight). In practice, when the N_eff is less than a certain predetermined threshold for the N_eff, some measures (resampling or redefining the covariance matrix in likelihood) are considered. Although the “proper” threshold depends on the problems and the purpose of the Bayesian updating, N/5 and N/3 are commonly used.

4. A case study

This section presents an example that applies PF to a typical geotechnical problem: deformation of the ground under embankment loading which is predicted by a numerical model. This example was used for a student competition on machine learning in 2020 organized by Technical Committee 304 and 309 of International Society of Soil Mechanics and Geotechnical Engineering (ISSMGE) [31]. Figure 9 shows the geometry of the cross-section of an embankment, which was constructed on layers of soft silty and clayey soils, installed with different types of instrumentation. In the figure, “M” and “VWP” indicate instruments to observe ground settlement and excess pore pressure of the ground. “HPG1” is an instrument to measure surface settlement. In this construction site, prefabricated vertical drains (PVDs) were installed to improve the ground because the ground consists of very soft clay layers. Information on the construction sequences of the embankment is shown in Figure 10. In the figure, four types of data measured at different locations (SP1 ∼ 4) were presented. In the student contest, the unit weight of embankment fill was assumed to be 20 km/m³. Deformation of the ground (settlement) was recorded at the site for around 3 years, which is shown in Figure 11. In the contest, only the settlement data measured for one year were presented, and the participants tried to predict the settlement after one year. The task of the contest was to develop a method to facilitate the prediction of the future settlement.

Figure 9.

Cross-section of embankment with location with instrumentation [31].

Figure 10.

Time history of embankment construction [31].

Figure 11.

Settlement data [31].

We used soil–water coupled Finite Element (FE) analysis with Cam-clay model for simulating deformation behavior of the ground. The FE mesh used for this problem is shown in Figure 12. We assumed horizontal displacements are zero at the edge of right and left of the model and vertical and horizontal displacements are zero at the bottom of the model. In addition, only the surface of the model ground is permeable: right, left, and bottom are assumed to be impermeable boundaries. Due to the symmetry of the embankment, only the right half was discretized for the numerical analysis. The foundation ground was divided into five layers; the first and fifth layers were assumed to be linear elastic bodies, and layers two to four were modeled as elastoplastic materials by the Cam-clay model.

Figure 12.

FE mesh used for the simulation.

In this study, a diagonal matrix with the maximum absolute value of the observation data used for PF is employed as the covariance matrix. This assumes that level of observation noise increases with increasing output (observation data). Such an empirical setting for the covariance matrix has been used in Bayesian updating with PF in geotechnical engineering [21, 22, 23], and despite its simplicity, there is a significant track record of obtaining reasonable results. The covariance matrix is given as follows:

$$R_t=\left[\begin{array}{cccc}1.08& 0& 0& 0\\ 0& 1.05& 0& 0\\ 0& 0& 0.74& 0\\ 0& 0& 0& 0.32\end{array}\right]$$E16

where R₁₁, R₂₂, R₃₃, and R₄₄ corresponds to the observation noise of M0, M1, M2, and M3, respectively.

In this case study, we focus only on the 2nd to 4th layers of clayey soil, which greatly influence the time-settlement behavior, and only the compression index λ and permeability coefficient k (m/d) of each layer are considered as the parameters for Bayesian updating. Therefore, the total number of geotechnical parameters to be estimated is six. Less influential parameters are assumed to be “deterministic,” and they are summarized in Table 2.

Layer	Soil type	Poisson’s ratio ν	Compression index, λ; Elastic modulus; E (Kpa)	Swelling index, κ	Void ratio e	Critical state parameter M	Permeability k (m/d)
1, 5	Sand	0.333	128,000	-	-	-	-
2,3,4	Clay	0.348	Unknown	0.123	1.525	1.023	Unknown

Table 2.

Soil parameters assumed to be deterministic.

While previous studies [21, 22, 23] mainly focused on the estimation of material parameters, this study also considers the embankment load and its loading history as unknown parameters and presents a specific example of load modeling. The reason embankment load and construction history have not been the focus previously is that nearly deterministic information can be obtained about the unit weight of the embankment and the construction process, which is considered negligible compared to the uncertainty of the foundation ground information. In this study, the uncertainty of the embankment load (f_t) at time t is modeled using the following equation:

$$f_t=\left\{\begin{array}{c}\mathrm{\alpha }{F}_0\left\{\left(exp\left(\mathrm{\beta }t\right)-1\right)/exp\left(\mathrm{\beta }T\right)-1\right\}\\ \mathrm{\alpha }{F}_0\end{array}\right.\begin{array}{c}t\le T\\ t>T\end{array}$$E17

where F₀ is the load at the completion of embankment construction, f_t is the load at time t, and α and β are parameters controlling loading history of embankment. Figure 13 represents the variation in load patterns due to differences in α when β is fixed at 1.0. Although the actual construction process shows a shape significantly different from an exponential function because of periods of rest in embankment construction, we proposed a new model based on an exponential function, which has relatively fewer parameters and is computationally easier to handle. Moreover, it is demonstrated through application examples that such a simplified modeling can function effectively in this case history. We added the parameters α and β into the unknown parameters, and the total number of parameters for Bayesian updating becomes eight. Since no prior information regarding the unknown parameters was provided, prior probabilities of the unknown parameters assumed to follow uniform distributions within the range shown in Table 3. During this process, several trial calculations were conducted to confirm the range within which the numerical results could explain the observational data, and then the parameter range was determined.

Figure 13.

Loading histories generated using Eq. (17).

	Parameter	Initial distribution
Loading model	α	0.05–0.50
	β	0.8–1.2
Cam-clay model	λ	0.4–1.2
	log10k(m/s)	0.01–0.50

Table 3.

Target parameters and their prior distributions.

Generally, it is known that the required number of particles N increases in proportion to the number of estimated parameters, but as mentioned earlier, the appropriate number of particles is often determined through trial and error. In this study, calculations were performed with 50,000 particles.

Figure 14 shows the impact of the presence or absence of load modeling on the Maximum A Posteriori (MAP) estimation of settlement. As evident from the figure, the load modeling enables the numerical analysis to successfully reproduce the observation data. On the other hand, without the embankment load modeled, a significant discrepancy between measurement and prediction is observed during the embankment loading period. This can be attributed to ‘modeling errors’ in the simulation model representing soil mechanical behavior. While altering the simulation model (Cam-clay model) may resolve this issue, the use of load modeling proposed in this paper leads to successful simulation of actual ground.

Figure 14.

Impact of the load modeling on settlement simulations.

Figure 15 shows the prior and posterior distributions for the eight unknown parameters. These histograms are drawn using 50,000 particles. Representing parameters as probability distributions allows for discussion of the uncertainty and reliability of the estimates, as well as the sensitivity of the parameters to the amount of settlement. For example, focusing on the compression index λ of the third layer, the posterior distribution is flat, showing little difference from the prior distribution. This indicates that this parameter has low sensitivity to settlement and that any value it takes does not significantly the result. On the other hand, the posterior distribution of the permeability coefficient k of the third layer is updated with a clear peak, indicating high sensitivity to settlement. In terms of sensitivity trends, the parameters of the load model (α and β) and the permeability coefficients k are highly sensitive, while the compression indices λ are less so. This can be attributed to the compression indices λ being in a trade-off relationship with the load intensity parameters (especially α) and due to the presence of multiple layers. The compression index λ is a parameter that defines the amount of settlement due to loading, and in the case of a single layer of clay ground, there are two possible combinations of parameters that can explain the amount of settlement. In other words, there are countless combinations of parameters that can explain the observational data, corresponding to the ‘ill-posed problem’ in the context of inverse problems, a challenge commonly encountered in analyses of inverse problems not only in geotechnical engineering but also in civil engineering.

Figure 15.

Prior and posterior distributions of unknown parameters.

Figure 16 shows a comparison between the settlement observed at M1 and the posterior distribution of settlement estimated by PF. These results consider load modeling, and hereafter, only results involved with load modeling are presented. The figure includes the 10%, 50%, and 90% percentile values. The percentile value represents the threshold below which 10% of the PF posterior distribution samples fall when ordered from smallest to largest. This allows for discussion of the accuracy of the prediction. Due to the settings in the simulations (covariance matrix Σ), the results of M0 and M1 have a significant impact, and the predictions match these two results. As a result, in M2 and M3, the measurement data do not fall within the 10% – 90% range, leading to an underestimation of settlement in M2 and an overestimation in M3. In this study, since the settlement at the ground surface is the most important for embankment structures, Bayesian updating focused on the prediction accuracy of surface settlement was conducted. However, if one wishes to focus on predictions for M2 or M3, adjustments can be made in the covariance matrix. By reducing the variance of the focused observation point, its influence becomes more significant in the Bayesian updating.

Figure 16.

Posterior distributions of settlement measured at different locations.

Figure 17 shows the posterior distribution of the settlement for M0–M3 as of June 1, 2016. The figure also displays each percentile value. The posterior distribution shows a nearly symmetrical Gaussian distribution shape, and irrespective of the measurement point, the 50% percentile value almost coincides with the mean of the distribution. It has been reported in previous literature that the posterior distribution being unimodal and symmetrical is common. This suggests the possibility of applying more efficient algorithms than PF in ground analysis, which is a subject for future research.

Figure 17.

Posterior distributions of settlement at June 1, 2016.

5. Conclusion

This chapter presented the theoretical and the practical effectiveness of the PF for geotechnical problems through a case history. First, some essential characteristics of geotechnical problems such as heterogeneity of the ground and variability of soil/rock parameters are explained to emphasize the motivation for using probabilistic methods in geotechnical engineering. Second, we have outlined the fundamentals PF by comparing two different sampling algorithms (SIR and SIS) and discussed which algorithm is more preferable in geotechnical problems from an aspect of geomechanics: stress history dependency of soils. Then, PF was applied to a real case history on embankment construction on clay ground to discuss the effectiveness in parameter (or probability density) identification and prediction in geotechnical engineering. The main conclusions of this study can be summarized as follows:

1. Probability estimation or uncertainty quantification plays a crucial role in geotechnical engineering because several different types of uncertainties (e.g., heterogeneity of the ground and variability of soil parameters) can be included in geotechnical design and construction process. Hence it is natural to use parameter identification methods that can provide posterior probability distributions of estimate such as KF and nonlinear KF in geotechnical engineering.

2. Unlike linear KF, the PF is applicable to nonlinear and non-Gaussian problems. Geotechnical problems usually show nonlinearity and non-Gaussianty, and the PF is preferable to apply in geotechnical engineering compared to standard KF.

3. Two types of sampling methods, SIR and SIS, are used in PF. Although SIR is commonly used in many different research fields, SIS is preferable to use for geotechnical problems because SIS can consider stress history soils undergo, which greatly impacts on mechanical behavior of soils.

4. PF can accurately estimate the posterior PDFs of geotechnical parameters in soil–water coupled FEM. In addition, PF can also be applied to identification of boundary condition, loading history in this study, for numerical simulations. The simulation with identified parameters predicts the future behavior of the ground with high accuracy: these results show that the PF has high potential for more general parameter identification and prediction problems in geotechnical engineering.