Parameter Identification Within Rocks Using Genetic Algorithms

This paper presents a concept for the application of particle swarm optimization in geotechnical engineering. For the calculation of deformations in soil or rock, numerical simulations based on continuum methods are widely used. The material behavior is modeled using constitutive relations that require sets of material parameters to be specified. We present an inverse parameter identification technique, based on statistical analyses and a particle swarm optimization algorithm, to be used in the calibration process of geomechanical models. Its application is demonstrated with typical examples from the fields of soil mechanics and engineering geology. The results for two different laboratory tests and a natural slope clearly show that particle swarms are an efficient and fast tool for finding improved parameter sets to represent the measured reference data.


INTRODUCTION
When compared to most other engineering tasks, geotechnical problems are often characterized by the following peculiarities.The materials involved are geological materials, that is, soils and rocks, which are inhomogeneous and consist of various phases in different states of aggregation.Initial and boundary conditions tend to be complex and heterogeneous.Furthermore, in real geotechnical field problems, the exact geometry is usually not known, with the available geometry information being limited to topographical surface data and punctual outcrops or soundings.For this reason, in geotechnics, there is always a need for a high level of simplification and abstraction.Frequently, continuum methods are used to calculate deformations in soil or rock, and the material behavior is simulated by means of constitutive models, which require a certain set of material parameters.
Normally, in geotechnical engineering, the values of these parameters are set based on the results of laboratory experiments, literature data, or even just experience values are used.The results of the calculation, a forward calculation, are then compared to measurement data obtained in the laboratory or in the field.Provided that the simplifica-tions made and the constitutive model chosen are appropriate and provided that the performed calculation gives plausible results, parameter values are then varied by trial and error in order to reach an improved fit of the calculation results to the measured data, the reference data.
Though this is done based on the experience of the geoscientist, the procedure remains to a certain extent arbitrary or at least subjective.
In recent years, due to the availability of sufficiently fast computer hardware, there has been a growing interest in the application of inverse parameter identification strategies and optimization algorithms to geotechnical modeling in order to make this procedure automated [1][2][3][4][5] and thus more traceable and objective.Furthermore, this approach provides statistical information, which can be used to quantify the calibration quality of the developed geotechnical model.
Applications of optimization procedures in geotechnics were described by many authors, for example, in the calibration process of geotechnical models [1,2], or to identify hydraulic parameters from field drainage tests [6].Already in 1996, Ledesma et al. [7] and Gens et al. [8] applied gradient methods to a synthetic and a real example of a tunnel drift simulation.Also during the excavation of a cavern in the Spanish Pyrenees the above mentioned group of authors was applying gradient methods for the identification of geotechnical parameters [9].Malecot et al. [10] used inverse parameter identification techniques for analyzing pressuremeter tests and finite-element simulations of excavation problems.For the identification of soil parameters, also genetic algorithms were studied [11,12].Feng et al. [13] used an inverse technique for the determination of the parameters of viscoelastic constitutive models for rocks, based on genetic programming and a particle swarm optimization algorithm.In the field of geoenvironmental engineering, Finsterle [14] examined the potential use of standard optimization algorithms for the solution of aquifer remediation problems in three-phase and three-component flow and transport simulations of contamination plumes.As a different aspect of parameter identification, Cui and Sheng [15] determined the minimum parametric distance to the limit state of a strip foundation by optimizing a reliability index.In 2006, J. Meier and T. Schanzin [5] applied particle swarm optimization techniques to geotechnical field projects and laboratory tests, namely, a multistage excavation and the desaturation of a sand column.
All cited references agree on the fact that back-calculation of model parameters by means of optimization routines is possible in the field of geotechnics, if an appropriate forward calculation depending on adequately realistic model assumptions is provided, for example, Calvello and Finno [1,2].In this context, particle swarms represent a powerful tool for finding parameter sets that best represent the reference data, with acceptable calculation effort and time consumption.

WORKING SCHEME OF THE ADOPTED PARAMETER IDENTIFICATION STRATEGY
The starting point of the parameter identification strategy presented in this study is given by an ordinary geotechnical modeling-task, the so-called forward calculation.This forward problem consists of a specified geometry with given initial and boundary conditions and a material model, which requires a set of material parameters to be determined.It is generally also possible to identify geometrical parameters [4], but this issue will not be discussed in this article.For the first run of the forward calculation, the user presets the unknown parameters, for example, entering estimated values, or the preset of the parameter vector is done by a random generator within values margins specified by the user.The relevant results of the forward calculation are then read out and their deviation from a set of reference data, usually measured data, is determined by means of an objective function.This procedure is repeated many times, while at any one time, an optimization algorithm, based on the parameter combinations and the values of the objective function during the previous forward calculations, identifies an improved parameter set to be used in the next forward calculation.This sequence of cycles, illustrated in Figure 1, is interrupted when one of the following stop criteria is fulfilled: (i) a maximum number of runs or maximum calculation time is reached; (ii) the deviation from the reference dataset, described by the value of the objective function, falls below a specified limit; (iii) the deviation could not be lowered during a certain number of cycles.
Hence we use a direct approach as described by Cividini et al. [16] to solve a back-analysis problem.In an iterative procedure, the trial values of the unknown parameters are corrected by minimizing an error function.It is therefore not necessary to formulate the inverse problem itself, the desired solution is obtained by combining the results of numerous forward calculations with an optimization routine.
To quantify the deviation between the reference data and the modeling results, we chose the frequently used and relatively simple method of least squares.In this method, the objective function f (x) for more than one reference dataset is defined as with (2) In ( 1) and ( 2), x denotes the parameter vector to be estimated and w g are positive weighting factors associated correspondingly with the error measure f g (x).Via the weights w g , the different series g can be scaled to the same value range and different precisions can be merged, for example, a series of measuring data possessing a higher precision is included with a higher weighting factor compared to more uncertain data.The particular numbers for the weights have to be given manually respecting the engineer's experience and they have to be specified depending on the optimization problem.The weighting factors w h are used to provide a possibility for considering different precisions and measurement errors within one and the same data series.The dimensions of the weighting factors can be taken in the way to obtain a dimensionless objective function quantity.For minimizing the objective function, we use a particle swarm optimization algorithm described by Eberhart and Kennedy [17,18].
A computer program developed by the first author of this paper implements this algorithm and disposes of interfaces to several commercial finite-element packages used in geotechnical engineering.A short pseudocode of the implemented PSO used in this study is presented in Algorithm 1.

CONCEPT FOR THE APPLICATION TO GEOTECHNICAL PROBLEMS
If the parameter identification strategy described above is to be applied to geotechnical tasks, the geotechnical model of the forward problem usually has to be adapted for its use in the optimization routine.The runtime of a single forward calculation has to be minimized in order to allow for a high number of calls.The number of calls needed depends furthermore on the number of parameters to be identified and, of course, on the used optimization algorithm.While the reduction of calculation time demands simplification and abstraction, the model still should be sufficiently complex to reproduce the reference data with the required accuracy.Furthermore, the number of required forward calculations can be reduced by applying hypersurface approximation methods [4].
As a next step, it is essential to select the parameters to be identified and to decide on the upper and lower limits of their plausible values margins.The values of some of the parameters might be fixed with the aid of previous knowledge.These specifications must be done with care and require the experience of the geoscientist, since they influence the obtained results.However, due to the application of a particle swarm optimizer instead of, for example, a gradient method, no initial guess for the parameter set is necessary because initial positions are generated randomly within the parameter value margins.
Due to the inhomogeneities of the geological materials involved and the uncertainties related to the initial conditions and the geometrical boundary conditions, geotechnical problems tend to be underdetermined.In order to improve and ensure the efficiency of the back-analysis, it is of significant importance to check if the set of parameters to be identified may be reduced and if for the prescribed trusted zone the optimization problem is well posed.For this purpose, a statistical analysis is done based on the results of a Monte Carlo procedure including a sufficient parameter set. Figure 2 shows the principle scheme of the matrix plot used here to visualize the results of the Monte Carlo simulation.A standard mathematical tool for examining multidimensional datasets is the scatter plot matrix (see [19]), whic is included in the matrix plot presented in Figure 2, where each nondiagonal element shows the scatter plots of the respective parameters.The matrix is symmetric.Matrix element D-B, for example, may suggest that the involved parameters B and D are not independent, but strongly correlated.The diagonal matrix plot elements (A-A, . . ., D-D) show plots where the value of the objective function is given over the parameter which is associated with the corresponding column.These plots are called hereafter objective function projections.If the problem is well posed, each of these plots of the objective function projections has to present one firm extreme value as it is the case in the diagram D-D.Otherwise, the respective parameter could not be identified reliably.By filtering out data points that have objective function values larger than a certain threshold level, the distribution of the remaining points gives a rough idea of the size and shape of the extreme value (solution) space.For further statistical analyses, the well-known linear 2D correlation coefficient can be calculated from the individual scatter plots.Referring to [20] in the analysis of our case studies, we consider variables with a correlation coefficient of less than 0.5 as "noncorrelated".

Description of the studied geological material and the adopted constitutive model
In the following examples of applications of the presented parameter identification technique using PSO, we are modeling the mechanical behavior of a natural soil.It is of geotechnical interest because it favors the development of numerous landslides, namely, rotational soil slips, earth slides, and earthflows.The studied material is clay that results from the weathering of structurally complex geological formations, the San Cassiano formation (Kassianer Schichten), and the La Valle formation (Wengener Schichten) of the Alpine Trias.These rock formations are made up by interbedded strata of marls, tuffites, claystones, limestones, dolomites, and sandstones.Like its source rocks, the soil is characterized by a high clay and silt content.In the field, it consists of a clayey matrix with coarser components, of diameters from centimeters up to meters, floating in it without mutual support.Therefore, it is not possible to sample the material as a whole representatively.As it has been completely remolded by earthflow phenomena, no preferred orientation of the components can be observed.For this study, only the fraction smaller than 2 mm was taken, as it is considered to determine the relevant mechanical properties of the entire soil and it  From this description of the material, it becomes clear that its mechanical behavior is expected to be very complex and therefore it is only possible to model some important aspects of this behavior.Like many soils with a high clay and silt content, the studied material is highly compressible and exhibits a significant amount of creep deformations, thus its behavior is strongly time-dependent.As a constitutive model, we chose the soft soil creep model, which was developed by Vermeer and Neher [21] to account particularly for these phenomena.The soft soil creep model requires the following material parameters to be specified (see Table 2).
A set of three parameters (c, ϕ and, ψ) is needed to model failure according to the Mohr-Coulomb criterion.Two further parameters are used to model the amount of elastic and plastic strains and their stress dependency.The modified compression index (λ * ) represents the slope of the normal consolidation line during one-dimensional or isotropic logarithmic compression.In the same manner, the modified swelling index (κ * ) is related to the unloading or swelling line.The modified creep index (μ * ) serves as a measure to simulate the development of volumetric creep deformations with the logarithm of time.
In the modeling examples of this article, we will not take into account the development and the influence of water pressures, which would be also possible but imply a considerable increase in calculation effort; and in the case of the slope example in Section 4.4, more reference data would be needed.All forward calculations were carried out applying the finite-element method using the commercial code PLAXIS (Version 8.2, professional, update-pack 8, build 1499) and considering the effect of large deformations by means of an updated Lagrangian formulation (updated mesh analysis) [22].

Oedometer test
A one-dimensional compression test was conducted by MFPA Weimar (Germany) [23] in a fixed oedometer ring with an inner diameter of around 7 cm (71.45 mm) and a height of around 2 cm (20.21 mm).Drainage was allowed on the top and at the bottom of the soil sample.All load steps were applied vertically, while the sample was held radially, impeding horizontal displacements.First, the sample was preloaded with 9 kPa during two days and with 13 kPa during one day.Then the load was doubled successively, with each load step lasting 24 hours, loading the sample with 25, 50, 100, 200, 400, and 800 kPa.After that, it was unloaded at 400, 200, 100, and 50 kPa, and finally it was reloaded again with 100, 200, 400, and 800 kPa (last step took 43 hours).The displacements of the sample top were recorded continuously.
For the numerical model of the test setting, we used an axisymmetric geometrical configuration with the exact dimensions of the test specimen.In order to minimize calculation runtime, the discretization was done with two six-node triangular elements only, which is the minimum possible number, as the software offers only triangular elements.Thereby, the duration of a forward calculation could be reduced to less than one minute on an ordinary personal computer.The accuracy of the deformation results was checked by carrying out comparative analyses with finer meshes.Horizontal fixities were assigned to the lateral boundary and to the rotation axis, simulating the stiff oedometer ring and vertical fixities were attributed to the basal boundary, representing the fix filter plate at the bottom.After generating the initial stress state by applying the soil self-weight (gravity loading procedure), distributed loads were applied perpendicular to the top boundary analog to the laboratory conditions.
Three parameters of the material model (λ * , κ * , and μ * ) can be determined directly from the oedometer test.As the material is known to show no dilatancy, ψ can be set to 0 • in all calculations.Laboratory data from shear tests on similar soil samples reported by Panizza et al. [24] is shown in Table 3.We averaged these values, giving double weight to the triaxial test data, which we assumed to be more precise, coming out with an average friction angle of 20 • and an average cohesion of 27 kPa.The set of experimental parameter values is shown by Table 4 and the results of a forward calculation using these parameters are presented in Figure 3 comparing them with the reference data of the oedometer test.
The graph shows that the deformations are underestimated by the simulation.In order to test the ability of the PSO algorithm to find good parameter combinations, wide search areas were chosen for the five parameters.A statistical analysis (see Section 3) comprising 2000 calls of the forward calculation was then performed varying these parameters.Their value margins are displayed in Table 5.
The scatter plot matrix of all data points with objective function values lower than 10 −7 is given in Figure 4.For the parameters λ * , κ * , and μ * , logarithmized margins of the search intervals were used in order to avoid overrepresentation of high parameter values.Furthermore, a parameter constraint was prescribed, demanding for λ * > κ * , which has to apply for all materials.The objective function projections of c, κ * , and λ * indicate that the data points showing good model fits seem to concentrate in quantifiable value ranges of these parameters, whereas μ * and ϕ cannot be identified.The modified compression index (λ * ) and the cohesion (c) appear to be correlated (correlation coefficient of 0.88), to a lesser extent, this holds true also for λ * and κ * (correlation   coefficient of 0.64).According to these findings, λ * , κ * , and c were selected for the optimization procedure.The friction angle (ϕ) and the modified creep index (μ * ) were fixed on their experimental values.After 159 cycles (1590 calls) the particle swarm optimizer had reduced the deviation to 5 * 10 −9 , which was found to be a sufficiently low value to stop the optimization routine.The identified best parameter set and the corresponding calculation results can be seen also in Table 4 and Figure 3.It becomes clear that the identified parameter set represents the measurement data much better than the available laboratory parameters.

Isotropic compression test
An isotropic compression test was performed in the triaxial apparatus on a cylindrical soil sample with a diameter of 5 cm (49.88 mm) and a height of 10 cm (98.55 mm).Drainage was allowed on top and at the bottom of the specimen.All load steps were performed isotropically, applying a hydrostatic cell pressure.The sample was preloaded at 30 kPa for three hours and, after that, at 50 kPa for 17 hours.Then it was gradually loaded to 800 kPa in one hour, increasing the load by steps of 100 kPa.After reaching this target load, the stress level was left constant for 3 weeks.Top displacements and volume change of the sample were recorded during the whole test.A reference dataset for the horizontal displacements was calculated from the vertical displacements and the volume change, assuming the shape of the specimen to remain exactly cylindrical until the end of the test.
For the numerical model of the test setting, an axisymmetric geometrical configuration with the exact dimensions of the test specimen was used.Again, the model was discretized only with two six-node triangular elements, to save calculation time.The accuracy of the deformation results was checked by carrying out several comparative analyses with finer meshes.Horizontal fixities were assigned to the rotation axis.Vertical fixities were attributed to the basal boundary, representing the fix filter plate at the bottom.After generating the initial stress state by applying the soil selfweight (gravity loading procedure), two independent and identical distributed loads were applied, one perpendicular to the upper boundary (vertically) and the other one perpendicular to the lateral boundary (radially).Loading was carried out the same way as in the laboratory, but instead of the stepwise application of the 800 kPa target load, this load was applied directly after the 50 kPa load step.For this reason, the 50 kPa load step in the model was prolonged in such a manner that the integral of the load as a function of time equals the test conditions.
In the example of the isotropic compression test, except for the relatively short phases before reaching the target load, only one load step (800 kPa) is applied.Therefore, of the model parameters only μ * can be determined directly from the test.This value (1.3 * 10 −3 ) is very similar to the one obtained from the oedometer test.Figure 5 shows the results of a forward calculation using this value together with the laboratory values of Section 4.2.
Also in this example, the deformations are underestimated by the simulation.A statistical analysis with 1820 calls was carried out varying the parameters c (cohesion), λ * (modified compression index), and μ * (modified creep index) within the boundaries given in Table 7, logarithmic values were used for the search intervals of the latter two parameters.
As the modeled test contains no unloading phases, it makes no sense to identify the modified swelling index κ * .Its value was therefore linked to the value of λ * by multiplying this parameter by 0.5, which is the typical κ * /λ * ratio, we observed in our laboratory tests performed on this material and similar materials.
The results of the statistical tests presented in Figure 6 suggest that good fits can be obtained for cohesion values between 20 kPa and 90 kPa, but apart from this, the cohesion value seems to have no influence on the quality of the model calibration.Whereas for the modified compression  index (λ * ) and the modified creep index (μ * ), the objective function projections suggest a preferred value range for the data points with low deviation values.Furthermore, it can be concluded that λ * and c might be quite closely related to each other (correlation coefficient 0.96), this means, for a given λ * , an appropriated cohesion value could be computed by the equation given also in Figure 6.This linear relationship seems to be valid for cohesion values between 20 and 90 kPa and λ * values between 0.064 and 0.165.
Therefore, only the parameters λ * and μ * were selected for the optimization procedure via PSO.We stopped this procedure after 500 calls (50 cycles).The identified values are shown in Table 6.In Figure 5, the calculation results are    compared to the reference dataset and the results obtained by using only laboratory data.
Like for the oedometer test, an improved parameter set could be found also for the isotropic compression test.It can be observed that the results are much better for the vertical displacements, although a weighting factor of one had been assigned to both datasets.
This may be due to the fact that the precision of the reference data is lower for the horizontal displacements, since the volume change of the sample could not be measured with the same accuracy as the top displacements.In addition to this, small inhomogeneities of the material or a small frictional resistance at the sample top could have caused a slight distortion of the cylindrical shape of the specimen which was assumed to calculate the horizontal displacements.

Analyzed section and reference data
The presented parameter identification technique was also applied to the 2D-model of a natural slope which is located in the municipality of Corvara in the Dolomites (Italy).It shows continuous creep deformations at the basis of a 20-40 m thick soil cover consisting mainly of the above studied material and very similar materials.As slopes of this type are known to show potential acceleration phases that can endanger human settlements, a detailed study of the geology and geomorphology as well as comprehensive monitoring were carried out by the University of Modena and Reggio Emilia and by the National Research Council's Group for Hydrogeological Catastrophes Defense (CNR-GNDCI) [25] by order of the Autonomous Province of Bolzano/South-Tyrol [24].In this context, the displacements of several surface points were observed regularly by means of global positioning system (GPS) measurements.In our study, a subarea of the slope was modeled along a representative 2D section.A site plan with the location of the section and the GPS points is given in Figure 7, whereas Figure 8 depicts a field survey of this section and shows the abstracted geometry model.

Geotechnical model
The geometry was determined using all the information available, that is, a core drilling near section-point B, the local geomorphology, refraction seismics, and direct current resistivity (DC-resistivity) [24].As exposed in Figure 9, the vertical profile of the slope was divided into three layers interpreting various inclinometer profiles reported by Corsini et al. [25]; the illustrated one is located near sectionpoint B. The uppermost layer, the soft soil cover, which is showing little internal deformations, was only considered in the form of its weight acting on the intermediate layer, the shear zone.Therefore, the displacement vectors of sectionpoints C and F are presumed to be equal to those of sectionpoints B and E, respectively.The shear zone is assumed to be a thin, soft, and highly plastic layer, exhibiting a pronounced time dependency in its mechanical behavior.The third layer is given by the underlying weathered bedrock, which is supposed to be stable.The earth pressure at the foot of the slope was assumed to be slightly lower than the earth pressure at rest; this means that there is no support obtained from the soil layer further downslope.A description of the assumptions made for the foot load is shown in Figure 10.
The actual main detachment zone is modeled as an open crack.No tensile forces are acting across it onto the downslope section of the sliding body, which is moving as a whole.Around section-point A, the soil body below the secondary shear zone is assumed to move at the same velocity as point B. The material properties of the secondary shear zone, which is not subject of this study, were fixed to comply with this criterion.At present, the displacement rates observed along the slope are more or less constant, being superposed only by seasonal variations attributed to fluctuations of the groundwater conditions which are not modeled in this study.Therefore, our geotechnical model features displacement rates remaining constant with time.

Numerical model
Figure 11 shows the characteristics of the numerical model of the studied slope.A plane strain geometrical configuration with the real dimensions of the slope was used.The model was discretized with 1070 triangular six-node elements.To save calculation time, the number of elements was reduced by modeling only the uppermost 20 m of the bedrock layer.The upper and the lower layers were meshed with the automatic meshing procedure of the software and using a very coarse setting.A linear elastic material model was assigned to them.Finally, one forward calculation took approximately three minutes on an ordinary personal computer.The intermediate layer was meshed manually by predefining geometry points in order to assure a sufficiently fine mesh and a suitable orientation of the triangles in order to go against an excessive distortion of the element shapes by the calculated deformations.For the same reason, the updated mesh option was only used in the last three years of the simulation (which are compared to the reference dataset).The detachment zones were modeled by means of interfaces on both sides of their geometry lines.The interfaces were modeled with a Mohr-Coulomb material model, with negligible cohesion, the same friction angle as the basal shear zone and with a constant reference stiffness in the order of magnitude of the shear zone stiffness.The accuracy of the calculated deformation results was checked by carrying out several comparative analyses with finer meshes and also with a horizontal basal boundary.Horizontal fixities were assigned to the lateral edges of the model, which extends over a total length of 1080 m.Vertical fixities were attributed to the basal boundary, representing the stable bedrock.

Calculation phases
In the first calculation phase, all three layers are made up by the bedrock material.An initial stress state is generated by applying the self-weight of this material (gravity loading procedure).In the second calculation phase, the two upper layers are replaced by the weaker material of the soil cover.The third calculation phase marks the starting point of the   slope instability in the model.The shear zone material with time-dependent mechanical behavior is inserted and the horizontal load at the foot is set.After that, the model is left creeping with unchanged boundary conditions during a period of 33 years.As the loading history of the shear zone material is unknown, this time period had to be chosen arbitrarily to reach constant displacement rates as they are presently observed along the natural slope.

Results of initial model using parameters derived from experiments
In a first trial forward calculation, for λ * , κ * , and μ * , the parameters calculated from the laboratory experiments were used as input values.As the deformations along the shear zone are known to persist since hundreds or thousands of years [26], the shear strength of this zone has decreased to a residual value that is characteristic for the soils originated by the weathering of the San Cassiano and La Valle beds outcropping in the whole slope area.Therefore, cohesion was assumed to be negligible (0.01 kPa) and a friction angle of 10 • was adopted, according to the average slope inclination observed in nearby areas which were formed since the Late Glacial by the studied processes (earth slides and earthflows) and covered by comparable soil covers [27].
The stiffness of the uppermost layer was set equal to the stiffness modulus observed in the oedometer test during unloading and reloading between the load steps 400 kPa and 800 kPa.For Poisson's ratio of this layer, we used 0.35, a value that is considered to be characteristic of clayey soils.
The experimental parameter values are shown in Table 8 and the deformations calculated on their basis for the last three years of the creep phase are presented by Figure 12.The latter are only in the range of millimeters, and thus not representing the actual situation in the field, where between September 2001 and September 2004, displacements from several centimeters to several decimeters were measured.

Results of statistical analysis and optimization procedure
A statistical analysis was carried out (which will not be reported in detail here).One interesting finding of this analysis was that the friction angle and the modified creep index appeared to be closely correlated (coefficient of 0.92).The parameters λ * , κ * , and μ * ; the friction angle; and the stiffness of the uppermost layer (represented by its shear modulus G) were chosen for the optimization procedure during which they were varied within the intervals specified in Table 9.
After 82 cycles, each of them consisting of 10 forward calculations, the procedure was stopped because, from then on, the deviation could no longer be reduced significantly.The resulting parameter set is also given in Table 9. Figure 13 depicts the calculated deformations using the identified parameter combination, together with the displacement vectors of the GPS measurement points.
It can be observed that the identified parameter set is able to reproduce the field measurements qualitatively.Because of the simplifications made in the model, no exact fit of the displacement vectors is possible.The presented back analysis procedure gives one of a number of possible approximate solutions to the geotechnical problem and the result returned by the particle swarm optimizer can be seen as a parameter set that best represents the reference data.

CONCLUSIONS
A back analysis procedure for the identification of material parameters of constitutive models applied to geotechnical problems was presented.This procedure represents a direct approach based on the method of least squares, correlation analyses, and a particle swarm optimization algorithm.The applicability and suitability of the technique was demonstrated by means of three examples from the fields of soil mechanics and engineering geology.The studied material was a natural soil.Besides being way more objective and less arbitrary than the conventional trial and error procedure, the outlined method provides valuable information on the quality of the model calibration, the uniqueness of an obtained solution, or the determinateness of the problem.In all three examples, the particle swarm optimizer was able to identify an improved parameter set after a justifiable amount of forward calculations.Further research should also concentrate on the identification of the geometrical parameters of geotechnical problems.

Figure 1 :
Figure 1: Flowchart of the adopted iterative procedure.
experiment) Simulation with parameters from laboratory tests Simulation with parameters from optimization with PSO

Figure 6 :
Figure 6: Scatter plot matrix for the isotropic compression test.

Figure 7 :
Figure 7: Site map with the location of the section and the GPS points.

Figure 8 :Figure 9 :
Figure 8: Field survey along the section and derived geometry model.

Figure 10 :Figure 11 :Figure 12 :
Figure 10: Earth pressure assumptions made for the foot of the slope.

Figure 13 :
Figure 13: Calibrated slope model-comparison of modeled displacement vectors (blue) with reference data derived from the GPS measurements (brown).

Table 1 :
Properties of the studied soil.

Table 1 .
Unless otherwise expressly stated, all tests and classifications were carried out according to the German standard DIN.

Table 2 :
Parameters of the soft soil creep model.

Table 4 :
Parameters obtained from experiments and identified parameters using PSO.

Table 5 :
Search intervals for the parameters of the oedometer test.

Table 6 :
Parameters from laboratory tests and identified parameters using PSO.

Table 7 :
Isotropic compression test search intervals for the varied parameters.

Table 8 :
Laboratory values of the parameters used for the slope example.

Table 9 :
Search intervals and identified parameters for the slope example.