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
ELSEVIER

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 compu-
tation 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 tradi-

tionally 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 struc-
tures themselves. By using field instrumentation results to estimate geo-
technical 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 prob-
abilistic 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 back-
analysis problem becomes equivalent to the minimization of an objective
function which depends on the values of field observations. Classical func-
tion 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 mini-
mization 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 identifi-
cation procedure is able to provide information on the structure and cap-
abilities 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 excava-
tion 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 mea-

surements and model predictions are due to an error measurement. There-
fore 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 pro-

portional 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 var-
iance, 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 indi-

vidual 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 inclin-
ometer 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 mea-
surements 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 long-
itudinal displacement of measurement point n, v,, will be
._
n

v, =
c
i=l
Ei + A (11)

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
(12)
=dekbmn = 2 min(i,j)
m=l n=l

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
i=l
aili + B (15)

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

.;,41 = i: f: cov(a,, 4lrn4
m=l m=l n=l

i i min(i.i) (16)
= o2c 2 ImlnS*, = o2 c- 1;
m=l n=l m=l

If the integration constant is not fixed, the elements of the E, matrix are

(&)ijzmpl;+
(Y) (17)
m=l

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 pro-
blem. 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 mea-
surements. 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 inclin-

ometer, 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 measure-
ment covariance matrix should be used.

NUMERICAL IMPLEMENTATION

Minimization algorithm

The solution of the identification problem defined above requires the mini-
mization 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 par-
ticular 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 para-
meters [eqn (4)] is in general non-linear, the form of a minimization algo-
rithm 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 itera-
tion 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~ (28)

where Ax = (x* - x) and the matrix

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’ cor-
rection 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, ~1is also increased until a smaller value of the objective function is
obtained. There is not a common procedure for evaluation of p and some-
times 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 fol-
lowing 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
10 A. Ledesma et al.

this approach, the finite element method has been coupled to the identifica-
tion techniques described above.
We assume that the measured values are nodal displacements. The process
can be generalized to the case in which relative displacements between two
arbitrary points are considered.
The finite element method formulated in terms of displacements gives a
linear system of equations:

Kx = f (31)
where K is the global stiffness matrix, x the vector of nodal displacements,
and f the nodal forces vector. The K matrix can be computed as

K= B’DBdV (32)
J V

where B is the geometry matrix relating strains and nodal displacements and
D contains the constitutive law as a relationship between stresses and strains.
The vector of nodal forces is defined as

f=
J s
N’odS (33)

where N is the shape function matrix, and o, for an excavation problem, is

the vector of stresses acting on the excavation boundary, S.
The displacements can be computed from eqn (31) and then can be com-
pared with the values measured, in order to check the identification process.
Using the finite element method and the Gauss-Newton algorithm, the
procedure used to estimate parameters is as follows:

(4 Assume a set of initial parameters po.

(b) Calculate the stiffness matrix K.
cc> Solve the direct problem and compute displacements x, x = K-If.
(4 Calculate the differences between the measured and the computed
values, x-x*.
Evaluate a new set of parameters p by means of eqn (28).
Check convergence. If Ap is still large or the norm of (x - x*) has
not stabilized, the procedure starts again from step (b).
Computation of the sensitivity matrix

Note that computation of A = ax/ap is needed in step (e) - eqn (28). One
possibility for determining it is to apply a finite differences technique. Alter-
natively, an “exact” procedure to evaluate the sensitivity matrix, using the
 Estimation of parameters in geotechnical backanalysis - I 11

finite element approximation, is presented. It will be shown later that this

procedure is also more efficient in terms of computing time.
Deriving eqn (31) with respect to the parameters and rearranging, we
obtain:

It is assumed that the materials involved in the problem are linear isotropic
elastic. The parameters to be identified are the Young’s modulus E, the
Poisson’s ratio v and the coefficient of lateral pressure at rest Kc,. For the
Young’s modulus and the Poisson’s ratio [lo], af/ap = 0, and

(35)

since
8K
B@BdI’
dp=. J dp

Hence computing dK/dp is equivalent to find the stiffness matrix sub-

stituting the aD/ap matrix for D, which allows an easy implementation in a
finite element code.
For instance, in a plane strain problem:

dD 1 l-v V 0
-=
ilE (1 + v)(l - 2~) 0” l-v 0 (37a)
0 ?

2V(2 - V) 1+2.7? 0
8D

(
E
-= 1+2v2 2V(2 - V) 0 ) (37b)
au (l+G)(l-22y)2 0 0
-- (l-2v)
2

Note that if v is assumed known and is not being identified, aD/aE is

constant and has to be calculated only once.
For the identification of the K0 parameter, eqn (34) is used, resulting in:

_&g
8X
(38)
aK0 0

since for a linear elastic material, aK/aKo = 0.

12 A. Ledesma et al.

The excavation is simulated in one phase, by applying on the excavated

boundary the opposite nodal forces to the initial stresses. So, using eqn (33),
we obtain:

J
af
-
aK0
=
I
s
N’%dS
dK0
=
s
N’cr*dS
0 (3%

where erg = (&$, c$, 0), CJ~is the initial vertical stress at depth y and

0; = (o$O,O) (40)

Note that to calculate dx/aKo, the finite element routines used to find
nodal forces due to excavation are employed to compute af/G’& by chang-
ing only the initial stress field from CO to o$.
This procedure allows us the direct use of the finite element approximation
in the computation of the sensitivity matrix, which is expected to reduce the
numerical errors. Also it is generally more efficient than a finite differences
technique which requires to solve n + 1 direct problems to estimate y1deriva-
tives, whereas here no additional systems have to be solved if the model is
linear elastic.
It should be pointed out that the evaluation of A in step (e) of the iteration
procedure requires to K-’ to be known which has already been computed in (c).
Actually, steps (d) to (f) are very fast, especially if the Poisson’s ratio v is not
identified (that means, if it is assumed constant), and only the Young’s modulus
E, and the K. parameter are involved in the process. For this particular case,
K=Ek (41)
where the k matrix depends on the geometry and boundary conditions of the
problem, and it is constant in the iterative process. E is the Young’s modulus
which changes in each iteration. Using eqns (36) and (38), it can be seen:

ax x ax 1 1 af
FE=2 ’ -=Ek-
r3Ko -~Ko (42)

where k-l& is constant in each iteration.

RELIABILITY OF THE ESTIMATION

Useful information concerning the parameter identification problem is

obtained by factorization of the sensitivity matrix A by eigenvalues and
eigenvectors. This is due to the fact that the rectangular linear system
 Estimation of parameters in geotechnical backanalysis - I 13

Ax = AAp (43)

can represent a least squares problem defined by the objective function (8)
[24,25]. A general decomposition of a m x n matrix can be made in the form

A = UAV’ (44)

where A is a diagonal matrix containing the k - non-zero eigenvalues (also

called singular values) derived from the problem:

Avj = ,+j , A’ui = Xivi (no summation) (45)

where
Xi=Xj, i=j, i<k; Xi=xj=O, i>k, j>k

U and V are, respectively, m x k and n x k orthogonal matrices containing

the measurement-space and the parameter-space eigenvectors corresponding
to the non-zero eigenvalues. The “pseudoinverse” matrix of A (in the sense
defined above) can be computed using the procedure defined by Lanczos [24]
and Lawson and Hanson [25]:

H = VK’U’ (46)

and therefore,
Ap = HAx = VA-‘UtAx (47)

The eigenvectors can be used for a reparametrization of the problem as

follows:

Ap+ = V’Ap ; Ax+ = U’Ax (48)

Then eqn (47) becomes:

Ap+ = ,-‘A,+, i.e. Ap: = (1 /Xi) Ax+ (49)

Equation (49) shows how a linear combination of the measurements (given

by the U eigenvectors) is directly and uniquely related to a linear combina-
tion of parameters (given by the V eigenvectors) through the inverse of the
corresponding eigenvalues. Furthermore, the value of Ai will determine the
variance of the corresponding linear combination of parameters. If the mea-
surements are statistically independent, then the variances of the increment
of parameters can be directly evaluated:
14 A. Ledesma et al.

var(Apr) = ($)var(Axt) (50)

i

Therefore

var(Api) = $ ($)var(A$) (51)

J

This is, in fact, a lower bound, because of the linearity assumed in the cal-
culation of the variances. Using the decomposition of matrix A shown in eqn
(44), eqn (51) can be written as

varAp = (A’ A))’ (varAx+) (52)

where varAp and varAx are vectors of variances of Api and Axi, respectively.
The factorization of the sensitivity matrix A and the results shown above
can be generalized if the maximum likelihood criterion is used. In this case
matrix A’ instead of A is used, computed from the factorization:

AC,-‘A = A’Att (53)

and eqn (52) generalizes to

varAp = (AC-‘A)-‘(varAx+) (54)

Therefore, in a maximum likelihood approach, the covariance matrix of the

identified parameters C, can be computed as:

C, = (AC-‘A)-’ (55)

This matrix agrees with the inverse of the Fisher information matrix [26], and
represents the reliability of the parameters estimated in terms of their covar-
iances. It is again a lower bound of the variances, because of the linearity
assumed in its computation.
The information density matrix can be defined as Q = AH in the least
squares case or Q = A’H’ if the generalization to the maximum likelihood
criterion is used. The information density matrix represents the inter-
dependency of the observations established by the model and the redundancy
of the data. When the number of measurements coincides with the number of
parameters, Q becomes the identity matrix.
Matrices C, and Q are very useful to assess the reliability of the identifi-
cation result and to examine the relation between measurements and para-
meters defined by the structure of the model.
 Estimation of parameters in geotechnical backanalysis - I 15

SYNTHETIC EXAMPLE

Description of the problem

To illustrate the performance of the formulation presented above, a synthetic

problem based on a tunnel excavation case is presented. We suppose that the
material is linear elastic, homogeneous and isotropic, with a Poisson’s ratio
equal to 0.49 in order to simulate undrained conditions and a specific weight
of y = 20 KN/m3. The parameters to be identified are the Young’s modulus
and the K0 coefficient defined using total stresses. Figure 1 shows the finite
element mesh corresponding to this example, and the 12 nodes used as mea-
surement points. Horizontal displacements in points 1 to 7 could represent
values obtained by means of an inclinometer device, whereas vertical dis-
placements in points 8 to 12 could represent measurements from a extens-
ometer. Excavation is made in one step, and due to the symmetry, only half
of the geometry is considered.

I 135 nl I

Fig. 1. Finite element mesh of the synthetic example.

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 proce-
dure 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 Vertical movements

Point Displacement (cm) 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 Gauss-
Newton 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 hori-
zontal movement) and 12 (maximum vertical movement) is clearly shown,

IO0
0 0 J Gauss-Newton
103
t. D-----13 J Marquardt

p Marquardt
10’

10-l

1o-3 =

1o-5

lo-’

1O-9

lo-‘+, , , , , , ,
1 2 3 4 5 6 7
ITERATION
Fig. 3. Evolution of objective function in the iterative process using Gauss-Newton and
Marquardt algorithms.
18 A. Ledesma et al.

9
0.9

8 0.8
m
y”
%
w 0.7
7

m K, v=O
0.6
o---+ E pi=10
6
1 d
LL
“i
I/ c-----a K, pi=10
I.5

5 II
I

I I I I I I I

1 2 3 4 5 6 7 0
ITERATION

Fig. 4. Evolution of the parameters in the iterative process using Gauss-Newton and Mar-
quardt 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 mea-
surements 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

I MEASUREMENT 1 1
01

-0 1

MEASUREMENT 5
01

ex 0.3

2 02

= 0.1
iz
z 0

z -01
0
i= /]
9
g 0:
z 0.1
Y
I- 0
b -0’
,”
t
E 0.:

g 0.2

0.1

-0 (

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 deter-
mined, thus only two variables (E and Ke) are unknown.
(b) Twelve measurements are used, corresponding to the example related
above.
20 A. Ledesma et al.

Cc> Twenty four measurements are available, as shown in Fig. 6. Points

l- 15 correspond to horizontal displacements and points 16-24 refer
to vertical ones.
(4 Fifty five measurements are available: 24 of them are the same as
used in case (c), and the rest are horizontal and vertical displace-
ments 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) pro-
vides 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

10000

1000

P
100 ;
E
\
‘\ l .
\ b’
‘\
n- ----- KO 10
----
------ --a

I I I I
1
2 12 24 55
NUMBER OF MEASUREMENTS
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 loca-
tions close to the measurement point with a much higher level of inter-
dependence 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.
 El
22 A. Ledesma et al.

0A
. 212
24
55 measurements
measurements
measurements 7I

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 con-
sidered. 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 inclin-
ometer 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 o----o C, diagonal

5
1 2 3 4 5 6-
ITERATIONS

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 estima-
tion problems, can be considered as a particular case of this formulation.
The probabilistic framework requires the use of the measurements covar-
iance matrix, which has been computed for three instruments used in prac-
tical 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-Mar-
quardt 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 identifica-
tion procedure, which provides a consistent framework to solve this kind of
problem in a systematic manner, considering the error structure of the mea-
surement process. The effect of the number of measurements and the influ-
ence 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 devel-
opments 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.
26 A. Ledesma et al.

3. Maier, G. & Gioda, G., Optimization methods for parametric identification of

geotechnical systems. In Numerical Methods in Geomechanics (Edited by Mar-
tins, J. B.). Nato A.S.I. Series, Reidel, Boston, MA, 198 1.
4. Arai, K., Ohta, H. & Kojima, K., Estimation of soil parameters based on
monitored movement of subsoil under consolidation. Soils Foundations, 24(4)
(1984) 955108.
5. Gens, A., Ledesma, A. & Alonso, E. E., Maximum likelihood parameter and
variance estimation in geotechnical backanalysis. 5th Int. Conf. Appl. Statist.
and Probab. in Soil and Struct. Engng, Vol. 2, pp. 613-621, 1987.
6. Gioda, G. & Sakurai, S., Back analysis procedures for the interpretation of field
measurements in Geomechanics. Znt. J. Numer. Meth. Geomech., 11 (1987)
555-583.
7. Honjo, Y., Limanhadi, B. & Wen-Tsun, L., Prediction of single pile settlement
based on inverse analysis. Soils Foundations, 33(2) (1993) 126-144.
8. Eykhoff, P., System IdenttJ?cation. Parameter and State Estimation. Wiley, Chi-
Chester, 1974.
9. Cividini, A., Maier, G. & Nappi, A., Parameter estimation of a static geo-
technical model using a Bayes’ approach. Int. J. Rock Mech. Mining Sci. Geo-
mech. Abstr., 20 (1983) 215-226.
10. Gens, A., Ledesma, A. & Alonso, E. E., Back analysis using prior information.
Application to the staged excavation of a cavern. 6th Int. Conf. on Numer. Meth.
Geomech., Iconmig 88, 1988.
11. Arai, K., Ohta, H. & Kojima, K., Estimation of nonlinear constitutive para-
meters based on monitored movement of subsoil under consolidation. Soils
Foundations, 27(l) (1987) 35549.
12. Ledesma, A., Identification de parametros en Geotecnia. Aplicacibn a la exca-
vacion de tt’meles. Ph.D. Thesis, Technical University of Catalunya, Spain,
1987.
13. Akaike, M., A new look at statistical model identification. IEEE Trans. Auto-
mat. Control, AC-19 (1974), 716-722.
14. Kashyap, R.L., Bayesian comparison of dynamic models. IEEE Trans. Auto-
mat. Control, AC-22 (1977), 715-727.
15. Edwards, A.W.F., Likelihood. Cambridge University Press, Cambridge, U.K.,
1972.
16. Tarantola, A., Inverse Problem Theory. Elsevier, Amsterdam, 1987.
17. Carrera, J., State of the art of the inverse problem applied to the flow and solute
transport equations. Nato Res. Work. Adv. in Anal. Numer. Groundwater Flow
and Qual. Model, pp. 5499583. Reidel, Dordrecht, 1988.
18. Baram, Y. & Sandell, N.R., An information theoretic approach to dynamical
systems modelling and identification. IEEE Trans. Automat. Control, AC-23
(1978) 61-66.
19. Kovari, K. & Amstad. Ch., Fundamentals of deformation measurements. In
Int. Symp. on Field Meas. in Geomech., Zurich, Vol. 1, pp. 219-239, 1983.
20. Gioda, G. & Maier, G., Direct search solution of an inverse problem in elasto-
plasticity: identification of cohesion, friction angle and in situ stress by pressure
tunnel tests. Znt. J. Numer. Meth. Engng, 15 (1980) 1823-1848.
21. Fletcher, R., Practical Methods of Optimization. Vol. I, Unconstrained Optimi-
zation. Vol. 2, Constrained Optimization. Wiley, Chichester, 198 1.
22. Scales, L. E., Introduction to Non-linear Optimization. Springer, New York,
1986.
 Estimation ofparameters in geotechnical backanalysis - I 21

23. Marquardt, D. W., An algorithm for least-squares estimation of nonlinear

parameters. J. Sot. Zndust. Appl. Math., 11 (1963) 431-441.
24. Lanczos, C., Linear DtBrential Operators. Van Nostrand, London, 1961.
25. Lawson, Ch. L. & Hanson, R. J., Solving Least Squares Problems. Prentice-
Hall, Englewood Cliffs, NJ, 1974.
26. Bury, K. V., Statistical Models in Applied Science. Wiley, New York, 1975.