Dynamic equivalent soil characteristics identification using earthquake records

1.   Introduction

Previewing past earthquakes shows that site geological conditions has a considerable effect on devastations resulting from Earth seismic motions (Tavakoli and Kutanaei, 2015; Choobbasti and Kutanaei, 2017a; Kutanaei and Choobbasti, 2017). In 1985's Mexico city's earthquake, although the fault rupture of earthquake was located about 350 km away from the city, but a small part of the city (with an area of about 6 km²) which had a site with soft clay deposits, was devastated considerably (Rezaei et al., 2015; Choobbasti and Kutanaei, 2017b; Anvari et al., 2017). It has long been recognized that local site conditions and local effects are significant factors on earthquake ground motion on soft soil sites and the need of taking them into account becomes pronounced when an earthquake strikes a region (El-Hussain et al., 2014; Abbaszadeh et al., 2011; Kutanaei and Choobbasti, 2016a; Mashhadban et al., 2016a). In more recent years, the inclusion of site effects in code provisions has received a large attention from the engineering community (Chavez-Garcia and Faccioli, 2000). Different studies have been conducted in order to collect evidence of local site effects as well as the soil response (Sarokolayi et al., 2015, 2016; Kutanaei and Choobbasti, 2016b). In report from past earthquakes it is often mentioned that the building located on top of the hills and altitudes were faced more obvious damages compared with those located outside the hills. Among earthquakes in which the effect of topography has been observed we can mention Lambesc (1909), Ferioli (1967), Michoacan (1985) and Chili (1985) earthquakes (Choobbasti et al. 2015; Kutanaei and Choobbasti 2015a). Other studies like those by Nakamura (1989) and Lachet and Bard (1994) deal with characterization of the soil properties. Previous studies indicate that soil properties are the most important parameters for site effect analysis. These properties are the shear wave velocity (v_S), the damping (ξ), and the natural frequency (f_n) of the soil system. Consequently, we may conclude that in order to obtain reliable estimations of the soil response, and in particular of site effects, it is necessary to have a reliable method to estimate the values of shear wave velocity, and the natural frequency for the site under investigation. Techniques for soil property estimation can be grouped into in-situ (geotechnical and/or geophysical surveys) and laboratory (geotechnical and/or soil mechanics) methods. Values of the dynamic shear modulus for low levels of strain can be obtained both in the laboratory and in-situ. However, at present, soil damping can only be determined in the laboratory.

Earthquake record has been proved a useful tool for estimating in-situ characteristics of soils in previous researches (Ghayamghamian and Motosaka, 2001; Kurita and Matsui, 1998; Mikami et al., 2003). Loh et al. (1989) have estimated the possible effects of the source (corner frequency, high frequency decay rate of acceleration) as well as the local soil characteristics (depth of soil deposit, soil density and shear modulus) on the spectral characteristics of ground motions by developing a stochastic soil amplification function for a uniform soil layer and using data recorded by SMART-1 array on Lotung site in the Taiwan region. After that, a system identification (SI) method was applied to identify the hysteretic behaviour of soil deposits (Loh and Yeh, 1992). Glaser (1995) summarized a number of identification techniques used to estimate soil properties from strong motion recordings. Glaser and Baise (2000) have estimated dynamic properties (shear stiffness and damping ratios) of the soils at Lotung site from seismic vertical array measurements (input-output data sets) using both time-invariant and time-variant parametric modelling methods (system identification) and modelling the soil system by a simple lumped mass. Mikami and Sawada (2005) used the equivalent linear method for the identification of the soil properties during the 1995 Hyogoken-nanbu earthquake. The identification process was successfully conducted, and the stress-strain relationships of the soils at the liquefied site were obtained from different depths all at once. Pavlenko and Loh (2004) performed nonlinear identification of the soil response at Dahan downhole array site during the 1999 Chi-Chi (Jiji) earthquake. They found that reduction of the shear moduli in the soil layers did not exceed 5%, and the soil response was virtually linear. Mikami et al. (2003) utilized modified Kalman filter methods with local iteration for the identification of nonlinear and nonstationary soil characteristics at a liquefied site during the 1995 Hyogoken-nanbu earthquake. They used borehole-array strong motions data. An identification of nonlinear and nonstationary soil characteristics was performed successfully; and nonlinear restoring force-displacement relationships including progression with time were obtained. Harichane et al. (2005) presented an analytical, numerical (linear method) and experimental method for identifying soil profile properties by using system identification and two free field records. To demonstrate the numerical efficiency and the validity of this approach, two examples have been treated. They reported that this approach offers the capability for more complete and rigorous characterisation of sites serving as support for constructions at reduced cost compared with the classical approach using laboratory and in situ tests, when ground motion data from previous earthquakes are available. Assimaki and Steidl (2007) used genetic algorithm for the inversion of weak and strong motion downhole array data obtained by the Kik-Net Strong Motion Network during the M_W7.0 Sanriku-Minami earthquake. They reported that downhole array recordings, provide valuable information on the elastic and nonlinear in-situ material response under true seismic loading.

Nowadays, optimizations play an important role in many industrial procedures (Soleimani et al., 2011; Kutanaei et al., 2011; Tavakoli et al. 2014a; Choobbasti et al., 2014; Kutanaei and Choobbasti 2015b; Mashhadban et al. 2016b). Considering the increase in cost of industrial products along with shortage of pure material, the importance of optimization is now more pronounced. PSOA is a heuristic global optimization method and also an optimization algorithm, which is based on swarm intelligence. It comes from the research on the bird and fish flock movement behavior (Kutanaei and Choobbasti, 2013; Janalizadeh et al., 2013; Tavakoli et al., 2014b). The algorithm is widely used and rapidly developed for its easy implementation and few particles required to be tuned. The main idea of the principle of particle swarm optimization algorithm (PSOA) is presented; the advantages and the shortcomings are summarized (Kutanaei et al., 2012). Due to its many advantages including its simplicity and easy implementation, the algorithm can be used widely in the fields such as function optimization, the model classification, machine study, neutral network training, the signal procession, vague system control and automatic adaptation control. PSOA has very deep intelligent background and it is suitable for science computation and general engineering applications. Yuan et al. (2009) have demonstrated PSOA worked on seismic wavelet estimation and gravity anomalies as well. Song et al (2012) present the application of PSOA to interpret Rayleigh wave dispersion curves.

In this work, the equivalent soil profile characteristics are obtained using the combination of PSOA and site response method. Finally this approach is utilized to identify the soil profile characteristics of different sites by using strong ground motions data recorded during the Bam earthquake of December 26, 2003.

2.   Linear site response

The linear approach is the simplest approach to evaluate the ground response and against Kelvin's visco-elastic model using transformation duction. The time history of bedrock motion is converted to frequency domain using Fast Fourier Transforms (FFT). Transfer function is then applied to the each of frequency (Fourier series), to obtain Fourier series of surface/ground motion. Then, inverse FFT is applied to obtain represent surface motions in the time period (Seed and Idriss, 1982).

3.   Particle swarm optimization algorithm (PSOA)

As previously mentioned, the main idea of the PSO is to mimic the social behavior of birds, which are referred to as particles in the remainder. This is achieved by modeling the fight of each particle by using a velocity vector, which considers a contribution of the current velocity, as well other two parts accounting for the own knowledge of the particle and of the knowledge of the swarm about the search space. In this way, the velocity vector is used to update the position of each particle in the swarm.

PSOA was firstly proposed by Eberhart and Kennedy (1995) based on the population (swarm) of particles. Each particle is associated with velocity that indicates where the particle is traveling. If t be a time instant the new particle position is computed by adding the velocity vector to the current position

being x^p (t) particle p position, p = 1, …, S, at time instant t, ${v^p}\left( {t + 1} \right)$ new velocity (at time t+1) and S is population size. The velocity update equation is given by

for j = 1, …, n, where κ₁ is a weighting factor (inertial), κ₂ is the cognitive parameter and κ₃ is the social parameter that are set to 1, 1.5 and 2 respectively. ${\omega _{1j}}(t)$ and ${\omega _{2j}}(t)$ are random numbers drawn from the uniform distribution (0, 1), used for each dimension j = 1, …, n. ${y_j^p\left( t \right)}$ is particle p position with the best objective function (goal function) value so far and $\hat y\left( t \right)$ is a particle position with the best function value so far. PSOA can be described as follows:

1) Randomly initialize the swarm positions ${ X} =$ $\left\{ {{x^1}\left( 0 \right), \ldots ,{x^s}\left( 0 \right)} \right\}$ and velocities ${ V} = \left\{ {{v^1}\left( 0 \right), \ldots ,{v^s}\left( 0 \right)} \right\}$

2) Let t = 0 and ${y^p}\left( t \right) = {x^p}\left( t \right)$ , p = 1, …, s.

3) For all p in {1, …, s} do:

If $f\left( {{x^p}\left( t \right)} \right) < f\left( {{y^p}\left( t \right)} \right)$ then set ${y^p}\left( {t + 1} \right) = {x^p}\left( t \right)$ else set ${y^p}\left( {t + 1} \right) = {y^p}\left( t \right)$ .

4) For all p in {1, …, s} do:

5) Compute ${y^p}\left( {t + 1} \right)$ and ${x^p}\left( {t + 1} \right)$ , using equations (1) and (2). If the stopping criterion is true, then stop. Otherwise set t=t+1 goes to step 3.

4.   General procedure

An important goal for geotechnical engineers is the ability to estimate soil properties without the measurement process disturbing the soil mass. In as much as an earthquake can be considered non-destructive, the archetypal large displacement excitation is earthquake strong motion, and the resultant soil motions are routinely recorded. At sites with installed vertical arrays, both ground motions into the bottom of the soil layer of interest, and out of the top of this layer are recorded, as illustrated by the cartoon in Figure 1.

Figure  1.  Configuration of the system identification method

DownLoad: Full-Size Img PowerPoint

Given this known input propagating upward from depth, and the output at the top of each horizon of the soil column, the behavior of the soil at each interval between accelerometers can be modeled by inverse theory. If a suitable model is chosen to represent the system of interest, the estimated model parameters will correspond to important mechanical parameters of the chosen model system, such as stiffness and equivalent viscous damping ratio. The goal of system identification (SI) is to invert the recorded data to estimate a model of a system, providing the needed mechanical information of that system. In this case, the data are input and output pairs of recorded earthquake ground motions, and the system is the intervening soil layers. PSOA was used to system identification based on the least squares minimization technique. The error function (deference between surface earthquake record and bedrock earthquake record) is minimized according to the site parameters, in the frequency domain. Figure 2 shows the general procedure which has been used in the present work. The PSO algorithm uses the data generated by a linear approach.

Figure  2.  Flow chart of the general procedure in optimization

DownLoad: Full-Size Img PowerPoint

5.   Results and discussion

5.1   Verification

In order to test the numerical efficiency and the validity of the present SI method, we applied it to the identification of characteristics (thickness h, γ unite weight, damping ratio ξ and shear wave velocity v_S) of soil layer (Table 1) overlaying a rigid rock. In this case, the identification can be done with respect to the complex form or with respect to the moduli of the soil amplification function between free surface and bedrock. To apply this method, it must get an initial guess for the model parameters. These estimations as well as the identification results are given in Table 1.

Table  1.  Characteristics of uniform soil layer

Parameter	Actual	Identified
h (m)	10.0	9.8
ξ (%)	10	10
$v_{\rm S}$ (m/s)	200	198
γ (kN/m3)	18	18

 | Show Table

DownLoad: CSV

Figure 3 shows the impact of the swarm size and the number of iteration on the SI performance (convergence behavior) in term of root mean square error. As can be seen a swarm with 150 particles has acceptable performance. This figure also clearly indicates that the rate of decrease in root mean square error is rapid in the first iterations, whereas the changes were considerably low up to iteration number of 200. Hence there is no significant change in root mean square error after 200 iterations; this value was selected to be used for the number of iterations in simulations. The soil amplification functions corresponding to actual, identified parameters are compared in Figure 4. This figure shows an extremely good agreement between identified and actual amplification functions but the identified parameters (h, ξ, V_s) are slightly different to the actual ones. In other hand, the ratio h/V_s, converges exactly to the actual value 0.05, about 200 iterations.

Figure  3.  Effect of the number of iteration on the SI performance

DownLoad: Full-Size Img PowerPoint

Figure  4.  Amplification function of uniform soil layer

DownLoad: Full-Size Img PowerPoint

5.2   Application

5.2.1   Earthquake ground motion

The powerful earthquake of December 26, 2003 almost destroyed the city of Bam that is located in the southeastern part of Iran. The magnitude of the earthquake was M_W 6.6 (USGS), its epicenter was close to City of Bam in Kerman province, and the focal depth was estimated to be 8–10 km. The free field record obtained in Bam station (Figure 5) showed maximum PGA of 0.8 and 0.7 g for the east-west and north-south horizontal components, respectively, and 0.98 g for the vertical component.

Figure  5.  The acceleration time histories (longitudinal, transverse, and vertical components) recorded at the Bam station

DownLoad: Full-Size Img PowerPoint

5.2.3   Soil properties

Geotechnical drilling was performed near the Bam station, to determine the subsurface layering characteristics and then evaluate the effect of local site conditions on the strong ground motion recorded in the city. The soil profile in the Bam station contains sandy clay on the top, which overlays the dense sand and silty sand at the bottom. Unit weight of the soil varied from 16.8 to 21.2 kN/m³, and the maximum water content was about 9.5%. The values of SPT blow count was measured about 15 on the top 4 m of the site and increased up to higher than 50 with depth. It should be noted that the water table was not encountered up to the explored depth of 30 m in this site. The thickness of the alluvium at the Bam station is believed to be greater than 60 m. The shear wave velocity, v_S and soil profile, of the subsurface layers were measured from down-hole tests and shown in Figure 6. The value of v_S varied between 100 and 670 m/s and increased with depth, indicating high stiffness of soil formations at the study region.

Figure  6.  Shear wave velocity, v_S, and soil type versus depth (SC: Clayey sand, SW: Well-graded sand, SM: Silty sand, SP: poorly-graded sand, GP: Poorly-graded gravel, GM: Well-graded gravel)

DownLoad: Full-Size Img PowerPoint

5.2.4   Generation of ground motion records for bedrock

Instead of using real ground motion records, the real accelerograms were adjusted and scaled with available software of SeismoArtif to get the artificial accelerograms based on the target spectra given in 2800 code of Iran for the soil class I (Table 2). SeismoArtif is a software which is capable of generating artificial accelerograms matched to a specific target response spectrum using different calculation methods and varied assumptions. Synthetic accelerogram generation and adjustment method has been used (Halldorsson and Papageorgiou, 2005). One of the accelerograms is produced by scaling the ground motion record of Bam Station obtained during the Bam EQ. In addition, in Figure 7, the normalized response spectra of soil conditions based on 2800 code of Iran and normalized response spectra of generated artificial spectral accelerations are shown. The base and surface acceleration histories for bedrock and free field motions are presented in Figure 8.

Table  2.  2800 code of Iran for site classification

Site class	Description	Average shear wave velocities (m/s)
I	Rock	≥ 750
II	Very dense soil and soft rock	375–750
III	Stiff soil	175–375
IV	Soft clay soil	≤ 175 m/s

 | Show Table

DownLoad: CSV

Figure  7.  Comparison of normalized response spectra

DownLoad: Full-Size Img PowerPoint

Figure  8.  The comparison of bedrock and measured free field motion in Bam

DownLoad: Full-Size Img PowerPoint

5.2.5   Identification result

The bedrock and surface acceleration histories are used in the simultaneous identification of soil profile characteristics (layer thickness, damping ratio and shear wave velocity and unit weight) of sites by minimizing on the modulus of spectral ratios. The results of identification of soil profile characteristics of sites are presented in Table 3. The average shear wave velocities obtained from downhole test in BH1 and BH2 are 471 and 423 m/s, respectively. The average shear wave velocity obtained from system identification method is 502 and 412 m/s, respectively. This result reveals a good agreement between system identification method and downhole test. Moreover present SI method has very high convergence speed (30 s). The results indicated that the peak horizontal acceleration is about 0.55 g at the base compared to 0.78 g at the surface, i.e., an amplification factor equal to 1.4. These results showed that the motion at Bam site was amplified with respect to the base motion by about 40%. In addition, the effects of local site conditions can be seen in Bam site in term of earthquake damage. In the location shown in Figure 9, soil has amplified the earthquake ground motion and consequently increased the damage in north Bam.

Table  3.  Equivalent characteristics of soil layer

	h (m)	ξ (%)	vS(average) (m/s)	γ (kN/m3)
BH1	35	7	502	18.4
BH2	31	9	412	17.6

 | Show Table

DownLoad: CSV

Figure  9.  Damage comparison between (a) before Bam earthquake and (b) after Bam earthquake

DownLoad: Full-Size Img PowerPoint

6.   Conclusions

In the present paper, an approach using linear site response method and PSOA for determining soil profile characteristics is presented. We have illustrated the capabilities of an efficient PSOA for free field record inversion, employed to obtain the equivalent soil profile characteristics (thickness, damping ratio and shear wave velocity and unit weight). This procedure is utilized to identify the soil profile characteristics of different sites by using strong ground motions data recorded during the recent Bam earthquake of December 26, 2003. This approach offers the capability for more complete and rigorous characterization of sites serving as the support for constructions at reduced cost compared with the classical approach using laboratory and in situ tests, when the ground motion data from previous earthquakes are available. The results from identification also results in a better understanding of earthquake hazard.