Identification of parameters in nonlinear geotechnical models using extenden Kalman filter

Direct measurement of relevant system parameters often represents a problem due to different limitations. In geomechanics, measurement of geotechnical material constants which constitute a material model is usually a very difficult task even with modern test equipment. Back-analysis has proved to be a more efficient and more economic method for identifying material constants because it needs measurement data such as settlements, pore pressures, etc., which are directly measurable, as inputs. Among many model parameter identification methods, the Kalman filter method has been applied very effectively in recent years. In this paper, the extended Kalman filter – local iteration procedure incorporated with finite element analysis (FEA) software has been implemented. In order to prove the efficiency of the method, parameter identification has been performed for a nonlinear geotechnical model.


Introduction
The back analysis procedures require reliable and sufficient measurement data, a robust numerical model and an efficient method for the solution of the inverse (optimization) problem [1].However in situ measurements are susceptible to uncertainties due to instrument inaccuracies, environment and/or measurement noise or human handling.Furthermore, for geotechnical problems the back analysis employs direct approach based on iterative solution of a forward problem by means of a numerical approximation (e.g.finite element analysis) and therefore the exact model response is basically unknown.
The Kalman filter (KF) method was developed by R.E.Kalman in early 1960s [2].Since that time it has been widely employed in signal processing, mechanical systems, etc. as a very successful method for state estimation and parameter identification.In recent years, some applications of the KF method have also arisen in the field of geotechnical engineering for estimation of material parameters [3,4,6].
In this paper, the KF method for model parameter identification has been successfully implemented for applications in geotechnical engineering.The method is based on extended KF technique combined with FEA.The underlying idea relies on the parameter estimation using the extended Kalman filter (EKF) algorithm, a nonlinear version of the Kalman filter, which is described in detail in the subsequent section.The method has been successfully applied for different linear and nonlinear geotechnical problems.Observations necessary for the EKF estimation are provided through finite element analysis and they include simulated measurements required in geotechnical applications.A demanding task of incorporating the EKF algorithm within the FEA is solved through appropriate software coupling in order to provide a smooth bi-directional exchange of the required data both for the EKF execution a e-mail: tamara.nestorovic@rub.deas well as for the FEA.The novelty represents also the verification of the proposed method for model parameter identification by means of its performance with respect to different geotechnical problems, like retaining wall structure -the problem which has according to the authors' best knowledge not been treated from this point of view sufficiently detailed before.
A formulation of the Kalman filter along with finite element analysis requires two equations of a state-space model: the state equation and the observation equation.The estimation of parameters presented in this paper is aimed at identification of geotechnical model parameters which do not change in time.Thus, a stationary transition of state estimates is suitable for our purpose [3].
Stationary system is expressed by the following state equation contaminated with noise at the k-th observation data (imaginary time t k ): where x k is the state variable vector containing model parameters to be identified, dimension n × 1; w k is a system noise vector contaminating the state equation, dimension n × 1. System noise is treated as white noise having zero mean E(w k ) = 0 and a covariance matrix Most materials involved in geomechanical applications behave nonlinearly.Even for a finite element formulation, resulting in a discretized system of linear constitutive equations, the nonlinearities will arise when solving inverse problems.In general, the constitutive relations depend nonlinearly on the material and model parameters.The demand of our estimation task, the inverse formulation, is to find material parameters incorporated within the constitutive laws from e.g. the known displacements.The finite element formulation for a general static problem is given by a relation between applied loads and displacements: K is the stiffness matrix, dimension l × l; x is the state vector associated with the material model used, as in equation (1) x has dimension n × 1; u is the displacement vector, dimension l × 1; f is the load vector, dimension l × 1.
The model response is observed at a finite number of nodal points within the finite element model.Measurement data in our case can involve displacements, stresses, pore water pressures, etc.The choice of the type of measurements depends on the particular geotechnical problem and the available measurement tools.In this paper, displacements at a finite number of nodal positions are chosen to be compared with measurement data.The reason for this choice is that displacement (settlement) is the quantity of interest in many geotechnical problems and hence due to the well developed equipment for measurement of settlements it would be possible to obtain real data as input for the model parameter identification procedure.
In agreement with the measurement data, we take the similar quantities from the FEA as observation data.The observation vector at observed positions at time t k is related with the displacement vector in the following manner where: y k is the observation vector (measurement), dimension m× 1; T is the transfer matrix, dimension m × l; v k is measurement noise vector caused by uncertainties in measurement equipment, dimension m×1.Observation noise is also white noise having zero mean E(v k ) = 0 and a covariance matrix The stiffness matrix, K(x k ), is generally a non-linear function of the material parameters x k and for a compact presentation here we express the equation ( 4) in a more general form: (5) For simplicity, the system noise w k and the observation noise v k are assumed uncorrelated.Furthermore, it is assumed that the parameter uncertainty is independent.The same applies to observation uncertainties.Thus the covariance matrices Q k and R k are diagonal.
The state equation ( 1) and the observation equation ( 5) are the basic equations representing the state-space system.Our goal is to estimate the state vector of the nonlinear system based on the knowledge about the modelled system and observation availability which are both contaminated by uncertainties.If the system is linear, the standard linear Kalman filter is the best estimator to be used.For our case of highly nonlinear system (geomechanical structure), a nonlinear Kalman filter (extended Kalman filter) is required.Unlike the standard Kalman filter, the extended Kalman filter linearizes the system at the current state.

Extended Kalman filter (EKF) -local iteration procedure incorporated with finite element analysis
First we introduce definition of a priori estimate and a posteriori estimate.A priori estimate is the estimate of the state before observation data are available.Each time the data observation occurs, a priori estimate will be corrected to the new state, which is called a posteriori estimate, that better represents the system.In mathematical expressions, a priori state estimate (estimate) and its covariance matrix at time t k are estimated as follows given the observation data up to time t k−1 : Whenever observation data at time t k are available, a posteriori estimate and its covariance matrix can be calculated as follows: The estimation procedure should be initiated by best initial knowledge of the considered geomechanical model.As initial condition, the initial state estimate and the covariance of the estimation error should be available: The initial state, which consists of model parameters, is chosen based on engineering experience or preliminary examination of the structure.The closer the initial parameter values to the true ones are, the better it is to begin the EKF.But this is not necessary.The initial covariance matrix is assigned rather arbitrarily large due to the lack of confidence in the initial choice of the state vector.Before observation data of the model are available, the EKF propagates the mean and error covariance of the state through time.The time update equations of the mean and covariance of the state are calculated in the following way: x where xk is the mean of state, xk = E(x k ), and P k is the estimation error covariance matrix, The superscript '-' denotes a priori estimate, and the superscript '+' denotes a posteriori estimate.
As soon as the observation data calculated by FEA are available, observation update procedure is applied.We apply here the local iteration scheme for the observation updates in order to continuously get better a posteriori estimate without calculating time update.At time t k in the local iterative loop i of the EKF -local iteration procedure, the Kalman gain K k,i , observation update of the state x+ k,i+1 and estimation error covariance matrix P + k,i+1 are calculated according to the following equations: where H, which is termed sensitivity matrix, is composed of the derivatives of observation data with respect to the state variables.H cannot be calculated analytically in our problems because in the observation equation h(x) is not known.Instead, a numerical approximation of H will be performed using forward results from FEA.K g is the Kalman gain matrix.Derivation of Kalman gain formula based on the theory of least-squares estimation may be found in [5,7].y exp is the expected measurement.Expected data are measured field data for the geotechnical structure being studied.In this paper, measurement data are produced synthetically by finite element simulations assuming that the model parameters are already known.
In the method used here, the error covariance matrix of estimation P k is enlarged with a modification weight W in every global iteration to obtain fast convergence as proposed by Hoshiya et al. [3].The aim is to make the fictitious difference between the true state and the estimated state more significant as it can be seen from the function J in equation (8).J is named cost function, objective function or return function.
In a stable system, the EKF always attempts to minimize the cost function J.The state error covariance matrix P is enlarged at observation step k, which means that an amount of pseudo-error is added to the filter.This means we make EKF think that its current state estimate xk is too far away from the true state x k .The new Kalman gain K k , equation (7), will also be enlarged with enlarging P k .The updated observation equation in turn uses K k to scale up the innovation )) and adds it to the a priori estimate x− k .Thus, state variable estimate x changes more quickly, and this makes convergence rate of the filter faster.
An independent treatment of the FEA program provides a significant advantage to the proposed method, in the sense that any FEA program can be incorporated into the process after the communicative interface between the main algorithm flow and the FEA program is created.In this paper, the main algorithm flow and the communicative interface are realised via MATLAB computing language and Python scripts, whereas the ABAQUS FEA software is used for performing the forward simulations.
In order to begin the filter process, an initial value x 0 for the state vector along with the corresponding estimation error covariance matrix P 0 have to be assigned.Usually, the initial state values are taken to be the most likely parameter values that can be inferred from prior knowledge or preliminary geomechanical tests.Values of the corresponding estimation error covariance matrix explains how much confidence the filter designer may have on the geotechnical model under consideration.
The effects of the process noise w k and the observation noise v k are omitted from system equation and observation equation for the sake of simplicity.However, their corresponding covariance matrices Q and R remain in the equations of time update and observation update of the EKF process.These covariance matrices are written without subscript k, which indicates that they remain unchanged over iterations.
The calculation of observation vector derivatives with respect to the state values at the current state estimate, H k,i , The iterative procedure explained consists of a local iterative loop for calculating observation update and a global iterative loop for the filtering process.The local iterative loop iterates a certain amount of times according to user setting and the global iterative process will continue until a stop criterion is met.Possible stop criteria are: 1) cost function achieves a predefined minimum value; 2) absolute difference of two consecutive estimated state vectors is less than a predefined tolerance; 3) the specified number of iterations is exceeded.In this paper, the third stop criterion is selected since the number of iterations is usually small and it also serves the purpose of examining the filter process at some later time after convergence has been reached.

Identification of Mohr-Coulomb elasto-plastic constitutive parameters
Identification of Mohr-Coulomb constitutive model parameters is carried out employing an example of a 4 m high retaining wall structure as depicted in Figure 1.A 30 kPa pressure and a 500 kPa pressure are applied on the surface of the backfill over the lengths of 6 m and 3 m and the problem is considered as 2D plane strain problem.The backfill structure defined in this way is similar to the one discussed in [8], but with modified geometry, materials and loads.
The Mohr-Coulomb constitutive model of the soil backfill is assumed to have a friction angle φ = 37 • , a dilation angle ψ = 5 • and a cohesion c = 10 kPa.Before the soil undergoes plastic deformation it behaves linear-elastically with a Young's modulus of 50 MPa and a Poisson's ratio of 0.3.The concrete retaining wall is assumed to be linear elastic with a Young's modulus of 21.3 GPa and a Poisson's ratio of 0.2.Density of the soil is 1900 kg/m 3 and of the concrete is 3000 kg/m 3 .
The finite element simulation is performed by means of ABAQUS dynamic analysis module and it is carried out in two steps, each step lasts for 5 seconds.Pressure loads are applied instantaneously at the beginning of step 1. Gravity load and prescribed displacement u * are ramped over step 1 and step 2 as it can be seen in Figure 2.  Four-node linear plane-strain finite elements are used in this finite element model.Observation data are taken from vertical displacements at four positions that are marked with dark color dots in the soil region as shown in Figure 1.
Observation data is composed of the observation point vertical displacements at the end of step 2 when the response of the model has attained its steady state.Figure 3 presents the displacements obtained with the set of true constitutive parameters of the observed points over the simulation time.
It has to be pointed out that for the purpose of this example not all parameters involved in the Mohr-Coulomb constitutive model are subject to identification.The elastic parameters, Young's modulus and Poisson's ratio, and soil cohesion are assumed to be known, whereas the friction angle and the dilation angle, are parameters that are selected to be identified and therefore they form the state vector x = φ ψ , where friction angle φ and dilation angle ψ are assumed as known in order to generate the synthetic measurement data.Friction angle and dilation angle are 37 • and 5 • respectively.
The EKF estimation process is run with only one local iteration.One single setting for noise covariance matrices Q and R and weight factor W for both initial state cases has shown difficulty in obtaining fast convergence.Thus, these settings are adjusted according to each initial state case.
Converged parameter estimates after 6 iterations are shown in Table 1.

Conclusions
Good convergence of the parameter estimation process demonstrated on a geotechnical problem presented in this work has proved that the EKF incorporated with modern FEA program has prospects to be a reliable method for parameter identification in geotechnical engineering.The important advantageous feature of the EKF method, namely the utilisation of the measurement uncertainties, which cannot be avoided in large scale geotechnical structures, makes this method very attractive and powerful.In addition, EKFprocedure performs well even in the environment corrupted by noise and due to the underlying characteristic of the algorithm, the state-value vector estimate (i.e.identified parameters) represents the values that best minimize the estimation error covariance.

Fig. 3 .
Fig. 3. Displacements at observed points over simulation time

Table 1 .
EKF -converged parameter estimates after 6 iterations: initial case 2 Parameter Initial state Estimated state Deviation