The Passive Surface Wave Methods for Shallow Engineering Exploration Based on the ESPAC Technology

. ESPAC method is a rapidly emerging field of seismological research, which can reflect the physical properties of the Earth's medium. In the process of using the ESPAC method, sometimes the noise of the original data is relatively large, and the raw data of each seismometer needs to be preprocessed, including operations such as de-averaging, de-trending, re-sampling, normalization, and filtering. The selection of the normalized method and the selection of the bandwidth of the filter are particularly important, and it will produce the wrong result if not handled properly. This article attempts to use the extended spatial autocorrelation (ESPAC) method to extract Rayleigh-wave phase velocity dispersion curves from the vertical component of the seismic stations’ microtremors, and proposes feasible and effective solutions to the selection of the normalized method and bandwidth of bandpass filtering.


Introduction
The surface wave is an elastic wave discovered by British scholar Rayleigh in theoretical calculations in 1887 and later confirmed by observation. With the indepth study of surface wave characteristics, it is found that the dispersion characteristics of surface waves in layered media, that is, the variation of surface wave velocity with frequency, make surface wave exploration possible. The dispersion characteristics of surface waves in layered media are the basic principles of surface wave exploration. The sensitization of surface wave velocity to the elastic parameters of layered medium (longitudinal wave velocity, shear wave velocity, layer thickness, density), especially the sensitivity to the S-wave velocity [1], enables the surface wave to invert the S-wave velocity structure of the formation. The S-wave velocity structure is closely related to the bearing capacity of the foundation and the dynamic characteristics of the soil layer. It is the basic parameter for the seismic safety evaluation of the engineering site, so the surface wave exploration has great significance.
The S-wave velocity structure can reflect the physical properties of the Earth's medium and is an important basis for distinguishing the structural characteristics of the low velocity zone in the Earth's crust. And the frequency-wavenumber spectrum (FK) method and the spatial autocorrelation (SPAC) method are two dataprocessing techniques used to obtain S-wave velocity from microtremor-array measurements.
In the 1960s, the United States established a large seismic observation network (LASA) with a radius of more than 200 km. In order to extract seismic information caused by nuclear tests from observational data, Capon (1969) [2] developed a frequencywavenumber method (FK method), which was also used to study natural source seismic wavefileds.
The spatial autocorrelation (SPAC) method is the method that Aki [3,4] firstly extracted the passive source surface wave dispersion in 1957. However, based on the disadvantage of SPAC method, the circular-array configuration, Ling and Okada (1993) [5] and Okada (1994) [6] proposed the extended spatial autocorrelation (ESPAC) method, which makes it possible to design arbitrarily shaped arrays.
This article attempts to use ESPAC method to extract the Rayleigh wave phase velocity dispersion curves from the vertical component of the seismic station's microtremors.
The strong amplitude of a large earthquake existing in the seismic record is an interference for calculating the dispersion spectrum, and the normalization of the original data is to normalize the noise amplitude to eliminate strong amplitude interference. After testing, this paper chooses to normalize in the time domain, and this method is called running-absolute-mean normalization. This method will be illustrated in the next section.
At the same time, the original data also needs to be filtered, but it is usually unclear which frequency band should be selected. In order to solve this problem, spectrum analysis is suggested to performed on each seismometer to observe which band of the data amplitude spectrum has stronger energy. And the filter window can be selected as this band. Of course, if there is a large difference between the spectrum of one or several seismometers than other spectrums, the data of these seismometers can be considered as bad data. It is recommended to delete these data because the spectrum of each seismometer should theoretically be roughly the same.

Data Preprocessing
Seismic data preprocessing procedures include detrending, de-averaging, cutting, bandpass filtering, time domain normalization, etc.
The original seismic record sometimes has a certain tendency, and the trend needs to be removed. The same goes for de-averaging. The data processing flow is different in sequence, which will have different effects on the cross-correlation results. De-averaging and detrending should be done after cutting in order to get better data processing results.
Cutting is the process of cutting long time sequences into short time series to increase the number of superposition. Appropriate increase in the number of superposition can significantly improve the signal-tonoise ratio of the cross-correlation.
Normalization of the time domain is to normalize the noise amplitude to eliminate strong amplitude interference.
Bensen [7] introduced five normalization methods. The first and most aggressive method is called "one bit" normalization ( Fig. 1b), which retains only the sign of the raw signal by replacing all positive amplitudes with a 1 and all negative amplitudes with a −1. This method has been used in a number of seismic studies and proven to be very effective; The second method is the threshold clipping method, the threshold is equal to the rms of the noise recording amplitude. An example is shown in Fig.  1c; The third method, illustrated by Fig. 1d, involves automated event detection and removal in which 30min of the waveform are set to zero if the amplitude of the waveform is above a critical threshold. This threshold is arbitrary and its choice is made difficult by varying amplitudes at different stations; The fourth method is running-absolute-mean normalization (Fig. 1e), and it is widely used now. This method computes the running average of the absolute value of the waveform in a normalization time window of fixed length and weights the waveform at the center of the window by the inverse of this average. That is, given a discrete time-series , we compute the normalization weight for time point n as: (1) so that the normalized datum becomes . The width of the normalization window (2N + 1) determines how much amplitude information is retained. A onesample window (N = 0) is equivalent to one-bit normalization, while a very long window will approach a re-scaled original signal as N →∞. An example result of the application of this method is shown in Fig. 1f. Before using this method, it is necessary to select the length of the window with different lengths for testing. It is found that half the maximum period of the passband filter works well, but the other lengths can still achieve similar effects.
The selection of the filter band in bandpass filtering will be elaborated in the fourth section.

Extended Spatial Autocorrelation (ESPAC) method
In order to obtain the phase velocity from data sets of a micro-tremor-array measurement, we used the ESPAC method proposed by Ling and Okada (1993) and Okada (1994), which were developed from the SPAC method by Aki (1957Aki ( , 1964. The basic principle of SPAC: Given a set of passive source surface wave receiving points, one of which is at the center, the remaining points are equiangularly distributed on the circumference (Fig. 2), assuming the center point and the angular frequency received by any point on the circumference is . Assuming that the Fourier transform of the microtremor signal at point is , the auto power spectrum of this seismometer is: (2) As for point C and X, their surface wave signals are and respectively. Its spatial autocorrelation function is as below: (3) The spatial autocorrelation coefficient is defined as the spatial autocorrelation function on all sides average of the direction: (4) The integral result of the above formula can be expressed as: (5) Where is the zero-order Bessel function and is surface wave phase velocity. By fitting the azimuthally averaged spatial autocorrelation function to the Bessel function, the phase velocity can be obtained. In determining from equation (5) using the SPAC method, we must select one of the following conditions to obtain the solution; either is constant (condition 1) or is constant (condition 2). It is noted that condition 2 leads to better results than condition 1 because phase velocity is a function of frequency (Ling and Okada, 1993;Okada, 1994).
The ESPAC method is derived from the SPAC method by using condition 2 to remedy the constraint in the array design (Fig. 3), which allows irregular arrangement, and is found to satisfy the following equation in the least-squares sense, instead of equation (5): The following equation is used to determine the error in evaluating . Because there is only one unknown in this equation, it is easy to obtain a solution. In this article, we use the least-squares algorithm to obtain the results of error [8]: (7)

Example and Results
In order to illustrate the effectiveness of the bandwidth selection of bandpass filtering, this paper selects the seismic data of Tongxiang, China for processing. Before using the ESPAC method to calculate the spatial autocorrelation coefficients, the original data needs to be preprocessed, including bandpass filtering, runningabsolute-mean normalization, de-averaging, and detrending. This article uses the running-absolute-mean normalization method in order to eliminate the effects of extreme data. And The bandpass filtering is described in detail below.
The vertical component of the seismic stations' microtremors is segmented by 4096 data points. First, perform Fourier transform on each piece of data, and the spatial autocorrelation coefficients are calculated according to equation (4) and (6). The phase dispersion velocity curves of Rayleigh surface wave are obtained by fitting with the first-order zero-order Bessel function according to equation (7). Fig. 4 is the phase dispersion velocity spectrogram of Rayleigh surface wave without filtering.
From Fig. 4, the unfiltered spectrogram is disorganized and we cannot extract accurate Rayleigh wave dispersion curves. However, we have proposed a solution to the above problem-calculating the auto power spectrum of each seismometer by Equation (2). And the filtered frequency band range can be obtained from the distribution of the auto power spectrum. Fig. 5 is a spectrogram of the 13 seismometers.
From Fig. 5, it can be clearly seen that the data amplitude is stronger in the 1-10 Hz band, and the data in 1-10Hz is exactly what we need. Therefore, the bandpass filter window is chosen to be 0.1-15 Hz. Fig. 6 is the phase dispersion velocity spectrogram of Rayleigh surface wave with bandpass filtering at 0.1-15Hz.
From Fig. 6, by selecting the correct bandpass filter window for filtering, a clear phase dispersion velocity spectrogram of Rayleigh surface wave can be obtained, and it's a critical procedure for the acurate extracting of Rayleigh wave dispersion curve from Rayleigh wave phase dispersion velocity spectrogram.

Conclusions and discussion
In this paper, the array of 13 seismometers is used to process the Rayleigh wave dispersion curve based on the ESPAC method. In terms of data processing, the running average of the absolute value method is effective and is now widely used. From Figure 4 and 5, using the auto power spectrum of each seismometer to determine the range of the filter band can obtain a clear spectrogram, which is of great significance for the extraction of the Rayleigh wave dispersion curve.
Surface wave exploration can be divided into actice source method and passive source method according to whether artificial source excitation is needed. The main method used in active source method is MASW method (Multichannel Analysis of Surface Wave), and the main method used in passive source method is ESPAC method. Since the dispersion curve extracted by active source surface wave exploration is mainly concentrated in the high frequency band (shallow surface), passive source exploration can extract lower frequency information, so the joint exploration of the surface wave active source and passive source shows better superiority.
Therefore, in the future work, we will try to use a variety of surface wave methods to jointly explore methods, and hope to obtain better exploration accuracy and scope of application.