Estimation of Parameters in Geotechnical Backanalysis: I. Maximum Likelihood Approach

ELSEVIER Computers and Geotechnics Vol. 18, No. 1, pp. l-27, 1996 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0266-352X(95)ooo21-6 0266-352X/96/$15.00+ 0.00 Estimation of Parameters in Geotechnical Backanalysis - I. Maximum Likelihood Approach A. Ledesma, A. Gens & E. E. Alonso Geotechnical Engineering Department, ETS Ingenieros de Caminos, Technical University of Catalonia (UPC), c/ Gran Capitin s/n, 08034 Barcelona, Spain (Received 9 August 1994; revised version received 27 January 1995; accepted 30 January 1995) ABSTRACT The estimation of soil and rock parameters based on field instrumentation data is a common procedure in geomechanics. The use of system ident$cation and optimization techniques allows the performance of this type of analyses in a more rational and objective manner. In this paper a probabilistic formulation for the backanalysis problem is presented. The procedure described involves the evaluation of the measurement covariance matrices, which are derived for some geotechnical instruments used in jield instrumentation. The algorithm used to solve the mathematical problem of optimization is also presented, as well as its coupling to aJinite element code. The algorithm requires the computation of the sensitivity matrix, which can be evaluated “exactly” in terms of thefinite element method. Finally, a synthetic example, based on the excavation of a tunnel, is presented in which the elastic modulus E and the Ko parameter of the material are identljiedfrom measured displacements. The eflect of the number of measurements and their error structure is also discussed. INTRODUCTION Determination of parameters in geotechnical engineering has been traditionally carried out based on the results of laboratory or in situ tests [l]. There is merit, however, in using for the same purpose, observations and measurements carried out during the construction of the geotechnical structures themselves. By using field instrumentation results to estimate geotechnical parameters it is possible to take into account the large scale structure of the soil or rock which is outside the possibilities of the other procedures of parameter determination. 2 A. Ledesma et al. In fact, backanalysis of the behaviour of real structures has always been an important part of the best geotechnical engineering practice [2-71. In recent times, however, the adoption of system identification [8] approaches and optimization techniques has allowed a more systematic and rational approach to this problem. Although there have been significant contributions in this field using a deterministic point of view, there are important advantages in using a probabilistic approach [9,10]. In this paper, a probabilistic framework based on the concept of maximum likelihood, that represents a conceptual alternative to the more classical Bayesian approach, will be proposed. The aim of the approach is to compute the parameters that maximize the likelihood of obtaining the field measurements actually observed. The backanalysis problem becomes equivalent to the minimization of an objective function which depends on the values of field observations. Classical function minimization techniques can be applied. In this paper, the constitutive model that defines the behaviour of the material is linear. However, the resulting optimization problem is non-linear so that iterative function minimization procedures must be used. The backanalysis involving non-linear constitutive laws is outside the scope of this paper [I 1,121. It is important to note that the process of parameter estimation is carried out in the context of a specified model that includes geometry, boundary conditions and constitutive laws for the materials. Because of the complexity of real engineering situations it is almost unavoidable that a numerical approach is required in practice to define the model. The parameter identification procedure is able to provide information on the structure and capabilities of the model. The topic of model identification, however, is not treated in this contribution [13,14]. The paper presents the basic formulations of the backanalysis problem in a maximum likelihood framework. Appropriate covariance matrices for the field measurement errors for various geotechnical instruments are then computed. Because of the need of defining the geotechnical model numerically, the question of the coupling of the minimization procedure to the finite element method is addressed. Afterwards the reliability of the parameters obtained and its relationships with the measurement errors is discussed. Finally, the performance of the procedure is illustrated by means of a synthetic example. In a companion paper an application to a real case involving the excavation of a tunnel in overconsolidated materials is described. BASIC FORMULATION Assume that a deterministic model relates some unknown parameters, p, and a certain set of variables, x. The measurements are represented by x*. Estimation of parameters in geotechnical backanalysis - I 3 Then the differences between measurements and predictions of the model (x*-x) are considered as an error, that can be defined in a probabilistic manner. The best estimation of the parameters is then found by maximizing the likelihood, L, of a hypothesis, p, given some set of error measurements, (x* - x). The likelihood of a hypothesis [15, 161 is proportional to the conditional probability of x* given a set of parameters p: L = w-(x*/P) where k is a proportionality constant. This formulation has theoretical and conceptual advantages [ 171: - It is not necessary to define the probability of a hypothesis, which has become a controversial concept in probability theory. - It does not require the model to reproduce the true system exactly [181. - The model parameters are considered fixed, but uncertain due to lack of information. This allows to introduce in a systematic way the prior information on the parameters available. Examples of this may be seen in Cividini et al. [9] and Gens et al. [lo]. If the model is considered to be correct, differences between field measurements and model predictions are due to an error measurement. Therefore the probability of measuring x* given a set of parameters, p, is the probability of reproducing the error measurements x*-x. Assuming that probability distribution as multivariate Gaussian, it is possible to write: P(x* - x) = J&-J exp[ - i (x* - x)‘(C,)-‘(x* - x)] (2) where (x*-x) = the vector of differences between measured and computed values using a fixed model C, = the measurements covariance matrix, which represents the structure of the error measurements m = the number of measurements. ( )’ is used to indicate a transposed matrix, and the likelihood is now proportional to the value expressed by eqn (2). Maximizing L is equivalent to minimize the support function: S= -21nL (3) 4 A. Ledesma et al. Note that the likelihood is a function on the parameters, taking into account that a model relating state variables x and parameters p has been specified: x = WP) (4) It follows from the above remarks that S= (x*-x)Cxpl(x* -x)+ln(C,I+mln(27r) -2lnk (5) On the assumption that the error structure is fixed, only the first term of eqn (5) has to be minimized, and in this case the objective function becomes: J = (x* - x)‘C,-‘(x* - x) (6) The case in which the error structure is not fixed is outside the scope of this paper, and will be treated elsewhere. If the measurements are independent and its errors have a gaussian distribution of probability with the same variance, the matrix C, has the form: c, = a21 (7) where I is the identity matrix. In this case eqn (6) represents a least squares criterion and the objective function is reduced to the simplest expression: J = (x* - x)~(x* - x) (8) If the m measures are obtained from Y independent instruments with individual covariance matrices (C,), eqn (6) can be written as J=>:( XT - Xi)'(Cx),'(Xf - Xi) (9) i=l which is the most general form of the objective function that will be used in this paper. COVARIANCE MATRICES The objective function selected to identify parameters in a geotechnical framework involves the covariance matrix of the measurements, as shown in eqn (6). If an instrument performs independent measurements, then the covariance matrix will be diagonal, as indicated in eqn (7). Estimation of parameters in geotechnical backanalysis - I 5 However, in some instruments designed for linewise observations [19], the errors of the measurements are not independent. For instance, if an inclinometer device is used to measure horizontal displacements along a borehole, the value of the displacement j (and therefore its error) is based on all the previously measured displacements. In this section, the covariance matrices for three instruments performing linewise observations - sliding micrometer, inclinometer and deflectometer [19] - are obtained. In order to do that, it is useful to express the covariance matrix for each instrument as where 4 is a scale factor which represents the global variance of the measurements made using the instrument i, and (E& contains the error structure of the instrument which depends on the apparatus itself. Obviously, if the measurements are independent and have the same variance, E, is the identity matrix, as shown in eqn (7). Sliding micrometer It is assumed that the line over which the measurements are performed is divided in p sections of length li. The sliding micrometer will measure the change of length, ei, between adjacent measuring points. Then, the longitudinal displacement of measurement point n, v,, will be n ._ v, = c Ei + A (11) i=l where A is an integration constant. To compute the covariance matrix, it is logical to assume that the primary measurements Ei, are independent, with a variance u2. Taken the integration constant A as exactly known (we are interested now in the general structure of C,), an element of the covariance matrix can be written as (CJij = COV[Vi, Vi] = COV [ 2 Ei, 2 Ej] = 2 k COV(Ei, Ej) m=l ?I=1 m=l n=l =dekbmn = 2 min(i,j) m=l n=l (12) and the corresponding element of the matrix E,: 6 A. Ledesma et al. (Ex)ij = min(i,j) (13) If the variance of the integration constant A is taken into account, then (EJij = min(i, j) + (F) (14) Inclinometer The inclinometer measures slopes o’i in fixed points along a line, and it is used to compute the horizontal displacements along a vertical line. The value ai is assumed to be small. Then the displacement perpendicular to the line measurement can be computed as n un = c aili + B (15) i=l where Ii is the length between two consecutive points of measurement, and B is an integration constant that expresses the horizontal movement of the initial point. Assuming that the value of B is exactly known, the C, matrix for an inclinometer will be Cij = COV [Ui, Uj] = COV [ 2 Cimlm, k m=l .;, 41 = i: f: cov(a,, 4lrn4 m=l n=l i i min(i.i) (16) = o2 c 2 ImlnS*, = o2 c1; m=l n=l m=l If the integration constant is not fixed, the elements of the E, matrix are (&)ijzmpl;+ (Y) m=l (17) Deflectometer The deflectometer can measure slope changes Ki in some points along a line. Then the displacement perpendicular to this line is given by (18) Estimation of parameters in geotechnical backanalysis - I I where C and D are integration constants. Without these constants, the covariance matrix is in this case: and when the integration constants have their own variance, Invariance of J The covariance matrices described in previous section have been derived assuming that displacements are the variables used in the identification problem. An alternative approach is to choose directly as the basic measurement variables the values of Ei, ai and IFS. Whatever approach is used, the objective function to be minimized must remain the same. Suppose that there is a linear relationship between the two kinds of measurements. For instance, in the inclinometer case: Ui = Wijf3j (21) where W can be obtained from eqn (15). Assuming, for simplicity, B = 0: Wij = 2 .lmSmj (22) m=l The two objective functions available for the identification problem will be: JI = Ax’C,‘Ax ; J2 = Aa’(-$)Aa (23) where Ax = (x* - x) and Aa = ((x* - LX). Using eqn (21): J2 = Ax’($WWt)-‘Ax (24) 8 A. Ledesma et al. and it can be seen that (25) which coincides with the measurement covariance matrix for an inclinometer, as indicated in eqn (16). This proof can be generalized to any other instrument, taking into account the fact that c, = cov[u, ll] = cov[wa, wa] = WW’cov[a, a] = ww’(u21) = u2WW’ (26) Hence, Ji = J2 and the objective function is invariant with respect to using any linear combination of measurements. Usually, displacements are taken as state variables, in order to simplify the coupling of the identification method with a finite element code and, therefore, the appropriate measurement covariance matrix should be used. NUMERICAL IMPLEMENTATION Minimization algorithm The solution of the identification problem defined above requires the minimization of a suitable function. Without losing generality, we can seek the minimum of eqn (6) instead of eqn (9). This can be achieved by means of a wide range of unconstrained optimization algorithms, according to the particular expression of the objective function. In general it is necessary to choose between methods that only need evaluations of the function (i.e. downhill simplex used in [20]) and methods that also require computations of the derivative of that function. In general, algorithms using the derivative of the objective function are expected to be more powerful than those using only the values of the function. In this work, a Gauss-Newton algorithm has been adopted, because of its good convergence properties when seeking the minimum of functions such as eqn (6), and because the derivatives obtained are also useful in providing information on the reliability of the parameters identified. Assuming that the relationship between measured variables and parameters [eqn (4)] is in general non-linear, the form of a minimization algorithm will be defined in terms of an iterative procedure: P~+I = ok + AP/~ , J(P~+,) G J(P~ (27) Estimation of parameters in geotechnical backanalysis - I 9 where Apk is the advance vector in the parameters space, computed in iteration k using a suitable algorithm. The Gaus-Newton method [21, 221 is based on an expansion of the objective function into a Taylor series, and gives the value of Ap in each iteration according to: Ap = (Atc;'A)-'~tC;lA~ where Ax = (x* - x) and the matrix (28) is called the sensitivity matrix. If the number of measures is m and the number of parameters is n, the size of this matrix is m x n. In general Gauss-Newton algorithm tends to exhibit a rapid convergence in this kind of problem, although it can be unstable sometimes. In those cases, an improvement of the algorithm proposed by Levenberg and Marquardt [23] has been used, in which eqn (28) is changed into: Ap = (AT;*A+ $'A~C$AX where p is an arbitrary real number and I is the identity matrix. If p + 0, the Gauss-Newton procedure is obtained. As p increases, the parameters’ correction provided by eqn (30) becomes smaller and tends towards the gradient direction of the objective function. The algorithm has to decide on the value of p in each iteration, depending on the behaviour of the procedure. If the value of J becomes smaller, p decreases, reaching 0 (or a very small value) at the minimum. However, if J increases, ~1 is also increased until a smaller value of the objective function is obtained. There is not a common procedure for evaluation of p and sometimes trial identifications must be performed in advance. In the examples of this paper, an initial value pi = 10 was used and for the next values the following criteria was employed: pu,+ 1 = pu, * 10 if J,+ 1 > J,, otherwise P,+I = UO, wh ere n is the number of iterations. Coupling to the finite element method A general identification procedure for geomechanical problems should be associated with a numerical model relating measurements and parameters. In 16 A. Ledesma et al. For the analysis, the displacements corresponding to the parameters E = 10 MPa and K0 = 1 are assumed as input data. The values of these theoretical measurements are presented in Table 1. The identification procedure must be able to obtain those parameters from the information provided by the measurements. As it is a synthetic case, there is no measurement error and the measurements covariance matrix is assumed to be the identity TABLE 1 Computed Displacements Used as Measurements in the Synthetic Example Horizontal movements Point Displacement (cm) Vertical movements Point Displacement (cm) 1 0.316 8 -3.598 2 0.556 9 -3.905 3 2.038 IO -4.935 4 3.550 11 -6.321 5 4.550 12 -7.048 6 3.920 7 2.421 0 p=o 0 pi=10 Fig. 2. Contours of the objective function and paths of the iterative procedure. VI and V2 are the directions of the sensitivity matrix eigenvectors in the minimum. Estimation of parameters in geotechnical backanalysis - I 17 matrix. The effect of the error structure on the identification process will be considered later. Starting with the p. values: E = 5 MPa and & = 0.5 the identification procedure has achieved the minimum in 4 or 5 iterations. Figure 2 shows the iterative process on the objective function J(E,&), using both, the GaussNewton algorithm (p = 0) and the Levenberg-Marquardt algorithm (initial value of p = 10). Figures 3 and 4 show the evolution of the iteration procedure. It can be seen that, in this simple case, there are not too many differences between the algorithms. In general, the Levenberg-Marquardt algorithm is slower but more robust than the Gauss-Newton, especially if the number of parameters to be identified is large. Finally, Fig. 5 shows the coefficients of the information density matrix of some measurements. The importance of the measurements 5 (maximum horizontal movement) and 12 (maximum vertical movement) is clearly shown, IO0 t. 0 0 D-----13 J Gauss-Newton J Marquardt p Marquardt lo-‘+ , , , , , , , 1 2 3 4 5 6 7 ITERATION 103 10’ 10-l 1o-3 = 1o-5 lo-’ 1O-9 Fig. 3. Evolution of objective function in the iterative process using Gauss-Newton and Marquardt algorithms. 18 A. Ledesma et al. 9 8 m % w 7 6 5 m K, v=O o---+ E pi=10 1 d LL I/ c-----a K, pi=10 “i I I I I I I I I I I 1 2 3 4 5 6 7 0 ITERATION 0.9 0.8 y” 0.7 0.6 I.5 Fig. 4. Evolution of the parameters in the iterative process using Gauss-Newton and Marquardt algorithms. which means they are the observations that provide more information in this particular case. The study of this matrix can be useful to check which measurements are more important, and at what points the field instrumentation should be installed. It must be noticed, however, that the components of this matrix depend on the value of the parameters. Therefore, it is necessary to have some suitable prior information about the parameters involved in the problem. Influence of the number of measurements If the number of measurements in this latest example is increased, we can expect an improvement in the identification results. Obviously that will depend on the kind of information added, i.e. displacements of nodes far from the tunnel will provide little additional information on the problem. In order to show the effect of the number of measurements involved in the results of the identification procedure, four cases based on the above example Estimation of parameters in geotechnical backanalysis - Z 19 01 0 -0 1 01 x e 0.3 2 02 = 0.1 iz z 0 z -01 0 i= 9 g 0: z 0.1 Y I0 b -0’ ,” t E 0.: g 0.2 0.1 0 -0 ( I MEASUREMENT 1 1 MEASUREMENT 5 /] I 1 I I , I I I I I I I 1234567 8 9 10 11 12 Horizontal Vertical MEASUREMENT POINT Fig. 5. Coefficients of the information density matrix for some measurements in the minimum of the objective function. have been considered. All of them use a diagonal measurements covariance matrix: (a) Only two measurements on the excavation boundary are available: vertical displacement of point 12 and horizontal displacement of point A (Fig. 1). In this case the problem becomes exactly determined, thus only two variables (E and Ke) are unknown. (b) Twelve measurements are used, corresponding to the example related above. 20 Cc> (4 A. Ledesma et al. Twenty four measurements are available, as shown in Fig. 6. Points l15 correspond to horizontal displacements and points 16-24 refer to vertical ones. Fifty five measurements are available: 24 of them are the same as used in case (c), and the rest are horizontal and vertical displacements from points distributed on the excavation boundary. The estimation of parameters is achieved in a few iterations in all cases. Some differences can be seen in the reliability of the solution for a fixed measurement error, represented by the parameters covariance matrix. As was expected (Fig. 7) the variances of the parameters identified, computed using eqn (55), are smaller when 55 measures are used [case (d)]. The case (a) provides a direct solution (equal number of measurements and unknowns), but the biggest uncertainty on the results. Fig. 6. Position of the measurement points in the synthetic examples (c) and (d). Estimation of parameters in geotechnical backanalysis - I 21 E ‘\ l \ ‘\ n- ----- KO ---- ------ --a I I I I 2 12 24 55 NUMBER OF MEASUREMENTS 10000 1000 P 100 ; \ . b’ 10 1 Fig. 7. Standard deviations of E and KO for different numbers of measurements used in the analysis. Values are proportional to the measurement’s standard deviation. The coefficients of the information density matrix for measurements A and B have been plotted in Fig. 8(a) and (b). Of course, if only two measures are used, the information is concentrated at points A and B. However, if 55 measurements are considered, the information is distributed along the locations close to the measurement point with a much higher level of interdependence between observations. lnclmometer Contour Extensometer Contour HORIZONTAL MEASUREMENTS VERTICAL MEASUREMENTS Fig. 8(a). Coefficients of the information density matrix in the minimum corresponding to point ‘A’, for different numbers of measurements. 22 A. Ledesma et al. El 0 A 0 . 2 12 24 55 measurements measurements measurements measurements 7 I 1 15 16 30 31 39 LO 55 IncLtnometer Contour Extensometer Contour HORIZONTAL MEASUREMENTS VERTICAL MEASUREMENTS Fig. 8(b). Coefficients of the information density matrix in the minimum corresponding to point ‘B’, for different numbers of measurements. Influence of the covariance matrix The importance of using the complete measurements covariance matrix can be demonstrated by means of this synthetic example. To this end, only a set of 15 horizontal displacements from previous example (c) have been considered. They were perturbed with a noise which was generated taking into account the measurement process of the inclinometer device. It was assumed that the standard deviation of the slopes measured was ga = 0.005 rad, and that error generated by standard random techniques was propagated to the horizontal displacements. Two analyses have been carried out using this perturbed input data. The first one uses the full covariance matrix as defined in eqn (17). The second one considers E, = I. As the values used as measurements have the inclinometer error structure, the analysis performed using the full covariance matrix should give (on average) better results. The parameter estimation process should provide the values E = IOMPa, K0 = 1, and the results of the iterative process indicate that the complete covariance matrix provides better results in the parameters identified. In Fig. 9(a) and (b) the objective function using full and identity covariance matrices has been depicted. Note that the minimum obtained using the wrong matrix (E, = I) is far from the actual parameters point. When the full matrix is used, the error structure is well represented and the minimum obtained is close to the real one. Taking into account the formulation described in the Estimation of parameters in geotechnical backanalysis - I 23 Fig. 9(a). Contour lines of the objective function using the full measurements covariance matrix. 1 1.5 % Fig. 9(b). Contour lines of the objective function using the identity measurements covariance matrix. 24 A. Ledesma et al. above sections, this fact should be expected, although only in a probabilistic sense. Figure 10 shows the evolution of the iteration process, starting from E = 5 MPa and I& = 0.5, gives a result in 4 or 5 iterations. w8 7 6 5 o----o C, diagonal 1 2 3 4 5 6ITERATIONS 12 3 4 5 6 ITERATIONS Fig. 10. Evolution of the parameters during the iterative process using full and identity measurements covariance matrices. Estimation of parameters in geotechnical backanalysis - Z 25 In a real backanalysis problem, matrix C, can be evaluated using the information available from the measurement devices and the procedure described above. However, an alternative way to compute its value is to introduce it as a new parameter in the identification procedure. This method [lo] is outside the scope of this paper and is treated elsewhere. CONCLUDING REMARKS The problem of parameter estimation has been formulated in a maximum likelihood framework. This approach has some advantages, i.e. allows the introduction of the error structure of the measurements in a consistent way and provides an estimation of the reliability of the parameters identified. Furthermore, the least squares criterion, commonly used in simple estimation problems, can be considered as a particular case of this formulation. The probabilistic framework requires the use of the measurements covariance matrix, which has been computed for three instruments used in practical geotechnical problems: sliding micrometer, inclinometer and deflectometer. This approach results in a mathematical problem of minimization of a suitable objective function, that can be performed by means of any available optimization algorithm. In particular, Gauss-Newton and Levenberg-Marquardt algorithms have been used in this paper, as they have proved to be convenient in these cases. The identification process requires the computation of the sensitivity matrix A defined as derivatives of the measured variables (displacements) with respect to the parameters. A procedure to calculate that matrix using the finite element method in linear models is described. Finally, a synthetic example has illustrated the features of the identification procedure, which provides a consistent framework to solve this kind of problem in a systematic manner, considering the error structure of the measurement process. The effect of the number of measurements and the influence of the covariance matrix on the solution has been examined in this case. An example, involving the identification of parameters in a real tunnel excavation problem, is described in a companion paper. REFERENCES 1. Jamiolkowsky, M., Ladd, C.C., Germaine, J.T. & Lancellotta, R., New developments in field and laboratory testing. 11th Znt. Conf. SMFE, San Francisco, CA, Vol. 1, pp. 57-153, 1985. 2. Terzaghi, K. & Peck, R. B., Soil Mechanics in Engineering Practice. Wiley, New York, 1967.