Andrés José Alfaro Castillo

CIEES, Bogotá, Colombia. E-mail: alfaro@ciees.org

ABSTRACT

The assessment of local site effects is one of the most important subjects in Engineering Seismology. In order to perform an assessment, it is necessary to determine the S-wave velocity structure of the site. Additionally, in some basins, it is very important to know the deep sedimentary structure, due to the amplification phenomena of low frequency waves. There are several techniques to achieve this purpose; probably the most inexpensive technique is using the vertical component of microtremors measured with an array of seismographs. The phase velocity of Rayleigh waves is inverted to an S-wave velocity (Vs) profile using optimization techniques. Most of the time, least square methods have been applied in the inversion. Recently, heuristic methods have also been used for the estimation of the S-wave velocity structure from microtremor.

In this study seven arrays of microtremors in the city of Tsukuba city were performed, located to the NE edge of Kanto Basin, in order to estimate the deep S-wave velocity structure. The spatial autocorrelation method SPAC was used to determine phase velocity dispersion curves in the frequency range from 0.3-2.5 Hz. The determination of Vs profiles reached a depth of 750 m. Two methods were used to estimate the Swave velocity structure: Inversion method and a heuristic method via the combination of Downhill Simplex Algorithm with a Very Fast Simulated Annealing Method. Comparisons with Vs from the existent results from PS-logging tests at the center of the array showed the reliability of the heuristic method.

Key words: Heuristic Method, Annealing Method, Inversion, Kanto basin, Tsukuba, SPAC, shear wave velocity.

RESUMEN

En este estudio se desplegaron siete arreglos para microtemblores en la ciudad de Tsukuba (Japón), ubicada en la parte Nororiental de la cuenca de Kanto, con el fin de determinar la estructura profunda de velocidad de las ondas S. Para determinar las curvas de dispersión de velocidad de fase en el rango de frecuencias 0.3-2.5 Hz se utilizó el método de la autocorrelación espacial SPAC. La determinación de los perfiles de Vs alcanzó una profundidad de 750 m. Se utilizaron dos métodos para estimar la estructura de velocidad de las ondas S: un método de inversión y un método heurístico vía la combinación del método Downhill Simplex Algorithm con el método Very Fast Simulated Annealing. Las comparaciones de la estructura de velocidades Vs con los resultados existentes de pruebas PS de registros de pozo en el centro del arreglo demostraron la confiabilidad del método heurístico.

Palabras clave: Método Heurístico, Método de Annealing, Inversión, Cuenca de Kanto, Tsukuba, SPAC, Velocidad de Cizalla.

INTRODUCTION

The geometry of the subsoil structure, the soil types and the variation of their properties with depth, along with the lateral discontinuities and the surface topography can produce large amplifications of ground motion and increase the damage during destructive earthquakes. For this reason the accurate knowledge of the geometry and the Vs structure of alluvial–diluvial deposits and the basement are very important in microzonation studies.

The Vs structure is usually determined in the field by using conventional seismic prospecting techniques (reflection, refraction, boreholes) and in the laboratory through dynamic tests on soil samples. The use of conventional seismic exploration methods presents some difficulties when the deep sedimentary structure needs to be determined. For example, in reflection and refraction surveys, the use of artificial sources such as explosives or vibrators is necessary, practices that sometimes are not easy in urban areas. Furthermore, the dimensions of the required arrays are large according to the desired penetration depth and therefore it is difficult to be deployed in populated areas. Additionally, the cost of large scale deep geophysical prospecting is high, and for this reason, in most cases during site effect studies, the depth of the seismic basement is limited to a layer with Vs larger than 400 m/s (Engineering bedrock) and not the real very deep underground reflector of the incident waves (Vs > 3000 m/s, Seismic bedrock). Additionally, the cost for implementing a deep borehole is also high and the results are valid only for a single site.

Microtremors techniques have been accepted during the last decades as a really good tool for reconnaissance and research of both shallow and deep soil structures (Alfaro, 2005a, 2005b; Alfaro and Yokoi, 2005; Alfaro 2006). There are several techniques using microtremors, during the last years, however Horizontal Vertical Spectral Ratio techniques (HVSR) have been used all over the world after Nakamura’s classical paper (1989) to determine the soil’s predominant periods and dynamical classification soil (Alfaro et al. ,2001; Bhattarai, 2005 ); HVSR has prompted several discussions due to the lack of robustness in the theory (Horike et al., 2001), Arai and Tokimatsu (2000, 2004) however, developed a technique that allows the determination of Vs Structure by means of inversion of HVSR. They developed a complete formulation assembling surface waves to achieve the inversion. In this research, however a microtremor array technique was used, because it has a robust theory and it has been used all over the world, mainly in Japan.

The method uses the microtremor records obtained at stations deployed in a triangular array. The measurements are taken simultaneously at all stations, which are operating for a short duration of time. The analysis of the microtremor records is performed through the spatial autocorrelation coefficient method (SPAC method) introduced by Aki (1957, 1965) and established by Okada (2003). It is important to mention another technique for arrays analysis, Frequency-wavenumber spectrum method (F-K), developed by Capon (1969) and applied to microtremors by Horike (1985) and Okada (2003). The F-K method is used to estimate the dispersion curve of Rayleigh waves and the velocity structure. The weakness of this method is the need of simultaneous measurement with several stations.

The SPAC method is based on the theory of the Stationary Random Functions, according to which, microtremor is considered as a stationary stochastic process both temporally and spatially. In this study, microtremor measurements were performed at one site, representative from the geological point of view, in the city of Tsukuba, where information of a deep borehole of 1300 meters is available (Hayashi, 2005; Hayashi et al., 2005). The practical aim of this study is to estimate a Vs profile, especially for depths larger than 500 m, reaching bedrock depth. To examine the efficiency and the accuracy of the method the results are compared with the borehole Vs profile at the site.

THEORETICAL BACKGROUND

Conventional Spac Method

Aki (1957, 1965) presented a theoretical background for estimating phase velocities by means of the SPAC method. For this research Yokoi (2005a, 2005b, and 2005c) developed a set of computer codes following Okada (2003) and Morikawa et al. (2004) in order to obtain the S-wave velocity structure.

The phase velocity at frequency ω will be obtained as the argument of the Bessel function. In this study we used a software b_fit developed by Yokoi (2005a) to find the optimum value for phase velocity. That program uses the Levenberg-Marquardt method and some subroutines from Press et al. (2003). To perform the inversion two programs were used: disp_sma1 (Yokoi, 2005b) and surf96 (Herrmann and Ammon, 2004). disp_sma1 is a program to obtain the optimum underground velocity structure for a given dispersion curve of Rayleigh waves based on the Downhill simplex method combined with the simulated annealing approach.

Simulated Annealing Method

The Simulated Annealing method is based on the idea of thermodynamics in which a metal melt reaches to a low-energy state through gradual decrease in temperature (Metropolis et al., 1953). After that Kirkpartrick et al., (1983) applied the idea to optimization problems with an analogy between thermodynamics and optimization as shown in Table 1. The misfit to be minimized during inversion corresponds to energy in thermodynamics, and parameter changes due to variations of material state. This move of parameters is controlled by cooling schedule of the system with decrease temperature.

The algorithm of the Simulated Annealing Method is shown in Figure 1. First, a cooling schedule and an initial Temperature are defined, T0. Search areas for all the unknown parameters are also defined before calculation. Then, an initial model, m0, is randomly generated within the defined parameter spaces. The misfit, E (m0) for the initial model is calculated using equation (1).

Where C^c(t_j ) and C^o(t_j ) are calculated and initial phase velocities at a period of t_j. σ(t_j ) and N are the standard deviation and number of observations. Next, a random perturbation is added to the initial model to generate a neighbor solution, m1. The misfit order is calculated again, E(m₁ ) for the neighbor model. If the difference of the misfits of the two models, ΔE = [E(m₁ )-E(m₀ )], is negative, m1 becomes the present model. If the difference is positive (m1 is the worse model), m₁ is still chosen as the present model with a probability P = exp(- ΔE/T). Because of the temperature-dependent probability, a model with high misfit is frequently chosen at high temperature. At low temperature, a not suitable model is not often selected and only a good model becomes the present model. After these processes have been repeated for all the parameters at predetermined times, the temperature is decreased according to the cooling schedule. The present model can be modified to the neighbor model near the global solution by repeating the above processes. Although, the Simulated Annealing Method is one of various local search methods using a perturbation of the model parameters, the Simulated Annealing Method works as a global search method at high temperature, because it allows an incremental change on the misfit surface. However, works as a local search method at low temperature. This feature is different from pure random search methods, such as the Monte Carlo search method.

There are several algorithms used in the Simulated Annealing Method. In this paper the Very Fast Simulated Annealing ((VFSA) Ingber, 1989) is used. In the Very Fast Simulated Annealing, a perturbation is generated using equation (2).

Downhill Simplex Method

The Downhill simplex method developed by Nelder and Mead (1965) method, requires only function evaluations, not derivatives. According with Press et al. (2003) the method is not very efficient in terms of the number of function evaluations that it requires. However, the Downhill simplex method may frequently be the best method to use if the figure of merit is “get something working quickly” for a problem whose computational burden is small. In this study some subroutines from Press et al. (2003) were used.

OBSERVATION OF MICROTREMORS

The data of microtremors were recorded in the North-Eastern part of the Kanto Basin, in the city of Tsukuba (Japan), where the depth to bedrock was found to be about 600 m with a deep borehole. Observations were performed on July 27 and 28 2005; additional measurements were done on August 4 and 16 to improved some data sets of non acceptable quality. Observations included seven arrays: one with 29 m radius (largest side of 50 m); two with 115 m radius (largest side of 200 m); two with 290 m radius (largest side of 500 m) and two with 520 m radius (within largest side 900 m). The geometries of the arrays are shown in Figures 2 and 3. Hereafter, these configurations are called R and B arrays. The smallest arrays (50 m side) and the middle array (200 m side) are shown in Figure 2.

There were simultaneous observations at four sites, with vertical-component velocity type seismographs (VSE12-CC manufactured by Tokyo Sokushin Ltd.), with a natural period of To = 10 s, also the short period velocity seismometers (HS-1, manufactured by OYO-Geospace) with 2 Hz natural frequency. The microtremors were recorded by digital recorders McSeis-MT (manufactured by OYO Corporation) with a resolution of 24 bits for A/D conversion, finally converted with an analog band pass filter, which was from 0.1 s to 5.0s, the sampling frequency was 100 Hz. When applying the SPAC method it is only necessary to gather the vertical component to obtain Rayleigh waves and avoid any kind of interference due to the presence of Love waves, making the analysis easier (Aki, 1957). The records were synchronized with the time code generated by Global Position System (GPS) clocks. Three sets of data on 200 m array, 10 sets of data with 500 m array (total duration 2 hours) and 10 sets of data with 900 m array (total duration 2 hours) were gathered. To ensure to check the stability of environmental conditions and instruments, data was gathered for all sensors in a single place. After this, stability analysis sensors were deployed in the field, starting with the smallest array.

On the other hand, data quality depends on the amount of interferences from anthropogenic sources like vehicles; in this study 500 m arrays recorded the highest interference due to the location of B5 and R5 stations, close to an Avenue with intense traffic, particularly trucks. An example of waveforms can be observed in Figure 4, which corresponds to microtremors recorded with a 900 m array.

Using the SPAC method it is possible to assess SPAC coefficients, as shown in Figure 5, which are functions of distance and frequency. For low frequencies, SPAC coefficients have maximum values.

The next step in the analysis is to fit a zero-order Bessel function of the first kind for every frequency, this procedure implies the verification of the values that could be used for fitting. Figure 6 shows some examples of fitting Bessel functions. For certain frequencies it is possible to use data from various distances, in the case of low frequencies, however, it is possible to use only data from the largest array. The b_fit program by Yokoi (2005a) was used. By means of inversion it is possible to assess the dispersion curve; examples are shown in Figure 7.

DETERMINATION OF VS STRUCTURE

The phase velocity was inverted to Vs structure using two methods: The least squares method (Herrmann and Ammon, 2004) and the combination of Downhill Simplex Algorithm with Very Fast Simulated Annealing Method (Ingber, 1989; Yokoi, 2005c). The misfit to be minimized during inversion corresponds to energy in thermodynamics, and parameter change does to move of material state. This move of parameters is controlled by cooling schedule of the system with temperature decrease. Figure 7 shows dispersions curves obtained with sets of 50 m data array plus R200 m data arrays plus sets of R900 m array. Figure 8 shows one example of results using surf96 (Herrmann and Ammon, 2004), and Figure 9 shows the results obtained using: The Very Fast Simulated Annealing Method, which are in good agreement with existing PS logging and borehole data, and the results derived by Herrmann and Ammon (2004).. Results using both methods are suitable, however one advantage of the heuristic method via the combination of Downhill Simplex Algorithm with Very Fast Simulated Annealing Method is the possibility of assessing thickness and Vs, meanwhile in Herrmann and Ammon (2004) it is necessary to fix either the thickness or the Vs of each layer. Advantages of the surfer96 software include the calculation of: standard error (km/s); mean residual (km/s); average residual (km/s) and percent of signal power fit in percentage.

These results indicate that the combination of Downhill Simplex Algorithm with Very Fast Simulated Annealing Method is a promising tool in phase velocity inversion. Probably, this is also true for other inversion problems in Engineering Seismology.

DISCUSSION AND CONCLUSIONS

Among several geophysical prospecting methods, passive methods have the advantage they do not need either artificial sources that could disturb people or expensive drilling.

Judging from the waveforms obtained, mainly for station R5 and B5 in 500 m arrays, they recommend avoid large avenues or roads, because waveforms could include interferences and results could be not suitable. Also, it is advisable perform the measurements during hours with minimum interferences due to punctual loads, like heavy trucks, or follow the recommendation of Apostolidis et al. (2004) to locate the stations at least 50 m far from avenues or roads.

The most significant advantage, of estimation of Shear Wave Velocity structure, using arrays of long period microtremors, is that the method allows for the reliable determination of Vs profiles down to large depths (about 600 m) with relatively small apertures of the deployed arrays (900 m). This is significant for accurate soil response studies in large cities, where open free spaces suitable for deployment of large conventional arrays are difficult to find and high energy sources cannot be easily used. .

ACKNOWLEDGEMENTS

The author gratefully acknowledges the opportunity and the support given by the Japan International Cooperation Agency JICA, the International Institute of Seismology and Earthquake Engineering IISEE. Also thank to Dr. Toshiaki Yokoi, for his permanent support during the research, to Mr. Koichi Hayashi and Mr. Kunio Aoike from OYO Corporation whom cooperation during field measurement was priceless. Most of the figures were draw using gnuplot (Williams and Kelly, 2004).

REFERENCES

• Aki, K. (1957). Space and time spectra of stationary stochastic waves with special reference to microtremors, Bull. Earthq. Res. Inst. 35, 415-457.        [ Links ]

• Aki, K. (1965). A note on the use of microseisms in determining the shallow structures of the Earth's Crust. Geophysics. 30, 665-666.        [ Links ]

• Alfaro, A. (2005a). Estimation of the shear wave velocity structures using arrays of long period microtremors, Individual Studies by Participants at the International Institute of Seismology and Earthquake Engineering. Building Research Institute. Tsukuba, Japan, 41, 15- 28.        [ Links ]

• Alfaro A. (2005b). Aplicación de los microtemblores en la Ingeniería Sísmica, XVI Jornadas estructurales de la Ingeniería de Colombia. Sociedad Colombiana de Ingenieros. Bogotá.        [ Links ]

• Alfaro, A. (2006). Determinación de la Estructura del Subsuelo utilizando arreglos de Microtemblores. Revista de Investigación. 6, no. 1, 1-6.        [ Links ]

• Alfaro A., and T. Yokoi (2005). Determinación de la estructura de velocidades de ondas de corte mediante arreglos de microtemblores de largo período. XIII Jornadas Geotécnicas de la Ingeniería Colombiana. Bogotá.        [ Links ]

• Alfaro, A., L.G. Pujades, X. Goula, T. Susagna, M. Navarro, J. Sánchez, and J. A. Canas. (2001). Preliminary map of soil's predominant periods in Barcelona using microtremors, Pure Appl. Geophys. 158, 2499-2511.        [ Links ]

• Apostolidis P., D. Raptakis, Z. Roumelioti and K. Pitilakis (2004). Determination of s-wave velocity structure using microtremors and SPAC method applied in Thessaloniki (Greece), Soil Dyn.Earthq. Eng. 24, 49-67.        [ Links ]

• Arai H. and K. Tokimatsu (2000). Effects of Rayleigh and Love waves on microtremor H/V spectra, Proc. 12th World Conf. Earthq. Eng, paper 2232, CD-ROM.        [ Links ]

• Arai H. and K. Tokimatsu (2004). S-Wave velocity profiling by inversion of microtremor H/V spectrum, Bull. Seism. Soc. Am., 94, no. 1, 53-63.        [ Links ]

• Bhattarai, M. (2005). Seismic microzonation using H/V spectral ratios with single station microtremor survey, Individual Studies by Participants at the International Institute of Seismology and Earthquake Engineering. Building Research Institute. Tsukuba, Japan, 41.        [ Links ]

• Capon, J. (1969). High-resolution frequency-wave number spectrum analysis, Procc. IEEE. 57, 1408- 1418.        [ Links ]

• Hayashi K. (2005). Active and passive surface waves, International Institute of Seismology and Earthquake Engineering. Lectures Notes. 45 pp.        [ Links ]

• Hayashi K. T., Inazaki and H. Suzuki (2005). Buried channel delineation using a passive surface wave method in urban area, International Institute of Seismology and Earthquake Engineering. Lectures Notes. 25 pp.        [ Links ]

• Herrmann R.B., and C. J. Ammon (2004). surf96 from computers programs in seismology, surface waves, receiver functions and crustal structure, version 3.30. Department of Earth and Atmospheric Sciences. Saint Louis University.        [ Links ]

• Horike, M. (1985). Inversion of phase velocity of longperiod microtremors to the s-wave velocity structure down to the basement in urbanized areas, J. Phys. Earth, 33, 59-96.        [ Links ]

• Horike, M., B. Zhao and H. Kawase (2001). Comparison of site response characteristics inferred from microtremors and earthquake shear waves, Bull. Seism. Soc. Am. 91, no. 6, 1526-1536.        [ Links ]

• Ingber L. (1989). Very fast simulated reannealing, Math. Comput. Modeling. 12, 967-973.        [ Links ]

• Kirkpatrick S., C. D. Gelatt, and M.P. Vecchi (1983). Optimization by simulated annealing, Science. 220, 671-680.        [ Links ]

• Metropolis N., A. Rosenbluth, M. Rosenbluth. A. Teller and E. Teller (1953). Equation of state calculations by fast computing machines, J. Chem.Phys. 21, 1087- 1092.        [ Links ]

• Morikawa H., S. Sawada, and J. Akamatsu (2004). A method to estimate phase velocities of Rayleigh waves using microseisms simultaneously observed at two sites, Bull. Seism. Soc. Am. 94, no. 3, 961-976.        [ Links ]

• Nakamura, Y. (1989). A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface, QR of RTRI 30, no. 1, February, 25-33.        [ Links ]

• Nelder, J.A., and R. Mead (1965). A simplex method for function minimization, Computer Journal. 7, 308-313.        [ Links ]

• Okada H. (2003). The microtremor survey method, Geophysical Monograph Series No. 12, Society of Exploration Geophysicists. Tulsa. USA. 127 pp.        [ Links ]

• Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B. P. Flannery (2003). Numerical recipes. The art of scientific computing. Code CDROM v 2.11, Cambridge University Press.        [ Links ]

• Williams T., and C. Kelley (2004). Internet Site for Gnuplot - an interactive plotting program, http://www.gnuplot.info        [ Links ]

• Yamanaka H. (2004). Application of heuristic search methods to phase velocity inversion in microtremor array exploration, Proc.13th World Conf. Earthq. Eng. Vancouver, B.C., Canada. August, Paper No. 1161.        [ Links ]

• Yokoi, T. (2005a). Combination of downhill simplex algorithm with very fast simulated annealing method - an effective cooling schedule for inversion of surface wave's dispersion curve. Proc. Fall Meeting of Seismological Society of Japan. S16-08010854.        [ Links ]

• Yokoi, T. (2005b). b_fit. Program for fitting the spatial auto-correlation coefficient determined from observed data to Bessel function and determine the dispersion curve. International Institute of Seismology and Earthquake Engineering.        [ Links ]

• Yokoi, T. (2005c). disp_sma1 Program to obtain the optimum underground velocity structure for the given dispersion curve of Rayleigh wave based on the down hill simplex method combined with the simulated annealing approach. International Institute of Seismology and Earthquake Engineering.        [ Links ]