Estimation of wave velocity and analysis of dispersion characteristics based on wavelet decomposition

1. Introduction

Wave propagation is characterized by the phenomenon of dispersion, which can be interpreted as a change in the waveform in space due to the different phase velocities of the components of the wave, resulting in a change in the energy distribution. The complexity of the crustal medium, stratigraphic boundaries, and the usual presence of rock pore structures containing fluids in the subsurface lead to complex attenuation and dispersion patterns of seismic waves [1]. Besides, strong ground motions generated by massive earthquakes can cause the Earth’s near-surface to soften and even lead to liquefaction in extreme cases. These changes may lead to a decrease in shear modulus and an increase in viscous damping [2]. The complex nonlinear characteristics of geotechnical materials under seismic action further exacerbate the attenuation and dispersion phenomena of seismic waves [3]. By studying the regularity of seismic wave dispersion, the information of underground medium carried in seismic waves can be extracted. Therefore, it is necessary to study the dispersion of seismic waves. Dispersion attribute studies of seismic waves require seismic data in different frequency bands, and hence time-frequency analysis of seismic data is required to obtain different frequency components [4]. Many methods can realize the spectral decomposition, such as Short Time Fourier Transform (STFT), Wavelet Transform (WT), S-transform, and Wigner-Ville distribution [5]-[6]. At present, the dispersion properties of seismic waves are usually studied without signal filtering during the time-frequency analysis process, and thus the attenuation and dispersion patterns of seismic waves are often difficult to be accurately identified due to the complexity of the stratigraphic structure. Based on the above, this paper proposes a wavelet decomposition-based wave velocity estimation and seismic wave filtering method, and establishes a numerical simulation example to verify the effectiveness. Then, the seismic records of the Tohoku Earthquake are used as the actual strong earthquake examples for filtering and dispersion analysis. Finally, it is concluded that the seismic waveform dispersion analysis method based on wavelet decomposition can effectively filter seismic waves and analyze the dispersion in a specific frequency band.

2. Wave velocity estimation and seismic wave filtering

2.2. Wavelet decomposition-based filtering method

In this paper, utilizing wavelet decomposition as a means of signal processing, a seismic wave filtering method based on the energy point of view is proposed (shown in Fig. 1).

Fig. 1Flowchart of the method

Before signal decomposition, it is necessary to determine a more appropriate frequency range, in which the signal energy needs to account for more than 90 % of the total energy. The determination of the frequency range needs to take into account the characteristic period of the seismic site, which needs to be set near the natural frequency. Due to the complexity of the stratigraphic structure, the seismic wave signals are usually very non-stationary. Therefore, focusing the frequency range of signal decomposition to the vicinity of the intrinsic frequency of the site can improve the accuracy of seismic wave dispersion analysis. The seismic waves are then decomposed into a total of $m$intrinsic mode functions (IMFs) in the selected frequency band using wavelet decomposition. The energy of each IMF is then calculated as follows:

2

$E_i={\sum }_{t=1}^N{f_i}^2\left(t\right),P_i=\frac{E_i}{\sum _{\mathrm{i}=1}^m{E}_i},$

where $f_i\left(t\right)$ denotes the Fourier transform of the $i$th $(i=\mathrm{1,2}···m)$ IMF, $t$ represents time, and $N$ denotes the total number of data points. Then the energy share of each IMF $P_i$obtained from the seismic wave decomposition can be calculated. And then the representativeness of the signal components relative to the original seismic wave can be determined. Finally, the signal components with energy share $P_i$less than 0.1 are filtered and the remaining IMFs are combined into a new seismic wave signal. Through the above filtering process, the low-energy clutter components in the frequency band can be eliminated.

3. Numerical modeling example

The feasibility and accuracy of the method are verified by establishing a numerical simulation model using SAP2000 software. Filtering and wave velocity calculations are carried out by using the data from the Tohoku Earthquake. A remote boundary homogeneous soil layer with a thickness of 381 m is modeled with the parameter settings shown in Table 1. The theoretical value of shear wave velocity $C_s$for this layer is calculated using the following:

where $G$represents the shear modulus and $\rho$ represents the density of the soil layer. The remote boundary homogeneous soil model has a spring support at the bottom of the model, which can simulate the input of seismic waves by inputting shear force and forward pressure (shown in Fig. 2). The distance between the left and right boundaries of the model is very far. When the seismic wave at the bottom has already traveled to the top, the seismic wave propagating to the left and right sides has not yet reached the boundary, which circumvents the influence of the reflection of seismic wave at the boundary of the two sides of the soil layer.

Fig. 2Numerical model of soil layers

Fig. 3Schematic diagram of the vs calculation

Firstly, the shear force $F_{bt}$ which will be applied to the spring mount of the model is calculated based on the seismic wave record shown in Eq. (4), where $K_{bt}$ is the spring stiffness coefficient on the viscoelastic boundary, $C_{bt}$ is the spring damping coefficient on the viscoelastic boundary, $v_0\left(t\right)$ is the incident wave displacement, ${\dot{v}}_0\left(t\right)$ is the incident wave velocity, and $C_b$ is the shear wave velocity. The calculation of coefficients $K_{bt}$and $C_{bt}$ is shown in Eq. (5), where ${\alpha }_t$ is the viscoelastic boundary correction factor, taken as 0.5, $C_s$ is the shear wave velocity of the soil layer, $\rho$ represents the density of the soil layer (shown in Table 1), and $R$ is the distance from the wave source to the artificial boundary:

Then the $F_{bt}$ (as the ground excitation) is input to the bottom boundary of the soil layer for simulation, and the time curve data of point A (mid-point at the bottom) and point B (mid-point at the top) is extracted after the simulation. Using the method of Sections 2.1 -2.2, the estimated shear wave velocity of this formation can be calculated and finally compared with the theoretical value $C_s$to verify the accuracy of the method.

Fig. 4Time history of Fbt

Fig. 5Numerical simulation results of point A and B

The seismic records from the KiK-net station in Japan for the Tohoku earthquake are selected for calculation. Taking the seismic wave record MYGH041103111446NS1 as an example, the shear force $F_{bt}$ applied to the spring at the bottom of the model by the simulated seismic wave input is first calculated (shown in Fig. 4). The seismic response curves at the bottom midpoint and top midpoint of this soil model are shown in Fig. 5.

From the aforementioned parameters of the soil model, the characteristic period $T_g$ of the soil layer can be calculated to be about 1.088 s, and the intrinsic frequency is about 0.919 Hz. The initial frequency band is [0.5, 2.0] Hz interval, and it is calculated that the energy of the signal in the band accounts for 99.92 % of the total energy, which meets the requirements. Then the displacement response of point B is decomposed by wavelet decomposition and the energy share of each IMF is calculated in the frequency band of [0.5, 2.0] Hz interval (shown in Fig. 6). From this, the IMF components with an energy share greater than 0.1 can be selected and the signal reconstruction is performed to obtain the new response of point B (shown in Fig. 7). It can be seen that the spurious components in the response time history of the simulation model are removed. The response of point A and filtered response of point B will be used as the input and output respectively to calculate the shear wave velocity $v_s$. The calculated estimation value of $v_s$ is 1.465×10³ m/s, which is close to the theoretical shear wave velocity $C_s$of 1.400×10³ m/s.

Fig. 6Energy proportion of IMFs

Fig. 7Filtered response of point B

4. Dispersion analysis of seismic waves

Still using the records MYGH04 of the Tohoku earthquake as an example for processing and analysis. Seismic wave data in the north-south and east-west directions is used for dispersion characteristics study, respectively. Firstly, considering the complexity of the stratigraphy and the fact that soil movement due to seismic action can have an impact on the propagation of seismic waves, it is necessary to pre-process records from surface station. In order to reduce the errors of the station's sensors under strong seismic effects, the first 15 s of seismic data are removed. Then the seismic wave data from surface are decomposed into multiple intrinsic mode functions (IMFs) based on wavelet decomposition. Using the estimation method of $v_s$ in Section 2.1, each IMF component can be processed to obtain the corresponding time-varying $v_s$. The relationship between shear wave velocity $v_s$ and instantaneous frequency is then obtained by using the Hilbert transform to obtain the frequency corresponding to each time point. The relationship between shear wave speed and instantaneous frequency is then obtained by using the Hilbert transform to obtain the frequency corresponding to each time point. However, the relationship curve directly calculated by the above process has many spikes and clutter components, so it needs to be further filtered to make the result smooth. Filtered results are shown in Figs. 8-9.

Fig. 8Wave velocity vs and frequency in east-west direction of IMFs

Fig. 9Wave velocity vs and frequency in north-south direction of IMFs

The $v_s$ and instantaneous frequency calculated from each IMF are then processed to obtain the average values based on empirical modal decomposition (EMD) method. In short, based on EMD, the $v_s$and instantaneous frequency calculated from each IMF can be decomposed into multiple signal components again separately. The energy share of the newly decomposed signal components is used as a weighting factor to calculate the weighted average value of $v_s$ and instantaneous frequency corresponding to each of the original IMFs, respectively. The results are shown in Fig. 10. Once the data from MYGH04 in the north-south and east-west directions have been calculated, the results in both directions can be combined to obtain dispersion results for the site (shown in Fig. 11(a)). A smoother result can be obtained after fitting the results of average value (shown in Fig. 11(b)). The mathematical expression used to fit the curve is as follows:

where $y$ represents the fitted value of$v_s$, and$x$ represents instantaneous frequency. Besides, $a=$ 424.5, $b=$ 2.727, $c=$ 8.027, $d=$ 1596. It can be seen that the estimated $v_s$at the site increases with increasing instantaneous frequency in the [0, 20] Hz band. And after the 20 Hz band, the $v_s$will grow in a slower trend and gradually stabilize.

Fig. 10Weighted average values of vs and frequency from IMFs

a) Results in east-west direction

b) Results in north-south direction

Fig. 11Dispersion analysis results of the site MYGH04

a) Average value of $v_s$ and frequency

b) Fit curve of $v_s$and frequency

5. Conclusions

In this paper, a wavelet decomposition-based wave velocity estimation and seismic wave filtering method is proposed, and the seismic records of the Tohoku Earthquake are used as the actual strong earthquake examples for dispersion analysis. Based on numerical simulation results and the analysis results of seismic record MYGHO4, the following conclusions are drawn.

1) This wave velocity estimation method can accurately calculate the $v_s$of site.

2) The wavelet decomposition-based seismic wave filtering method can effectively preprocess seismic wave data, filtering out low-energy clutter components through signal decomposition and reconstruction.

3) The filtering and signal smoothing methods based on the energy perspective can effectively overcome the difficulties caused by the complexity of soil layer structure in analysis, providing a potential new approach for frequency dispersion research.