Detection of Long-period Seismic Waves

The significant features of the far distant earthquakes can be efficiently identified and extracted through the signal processing and feature extraction techniques, despite that the seismic ground motion data are inherently non-stationary in both the amplitude and frequency. In this study, the selected techniques can be divided into the singular spectrum approach and the wavelet approach. In the singular spectrum approach, SSA is capable of decomposing the time series into a sum of informative components without prior knowledge, hence it gives a truly non-parametric tool for recognition and identification of a time series structure. In the wavelet approach, WPT is able to decompose the time series in both the time and frequency with an adaptive resolution; therefore, it

doi:10.6342/NTU201901242 40

provides an effectively diagnosing tool for the frequency content and the local information. In order to identify and extract the ground motion features from the recorded seismic waveforms, both moving window SSA and WPT were utilized to detect the existence of long-period seismic waves after the initial P-waves arrival since these seismic waves have shown the ability to bring a significant response to the wafer scanners while the vibration amplitude is large enough.

3.3.1 Energy Entropy and Major Components

The study in Section 3.2 has shown that the major frequency of the far distant earthquakes is quite low and the long-period seismic waves can force the wafer scanners in the high-tech fabs to shut down without any warning. To further investigate the major frequency varied with time, WPT was used to analyze the seismic waveforms of the three far distant earthquakes, i.e., the Indian Ocean earthquake, the Samar (Philippine) earthquake, and the Lushan (China) earthquake, recorded from the B550 station (from the NCREE network). Only UD direction was conducted because the P-waves is the fastest of seismic waves. In WPT, the Biorthogonal 6.8 was selected as the wavelet function, and 11 was assigned as the reference level of decomposition to provide high resolution in the frequency domain. The wavelet packet components can be reconstructed using Equation (2.33) and the amplitude function of these components can be estimated by Hilbert transform (HT). This transform derives an imaginary counterpart for a real-valued time series such that the resulted complex function is an analytic representation of the time series (Hilbert, 1953; Johansson, 1999).

𝒜(𝑡) = 𝑥(𝑡) + 𝑖HT[𝑥(𝑡)] (3.3)

where HT[∙] is HT operator that produces a harmonic conjugate of a real-valued time series, 𝑥(𝑡).

Theoretically, it is defined as an integral that can be evaluated by the Cauchy principal value:

HT[𝑥(𝑡)] =_𝜋¹∫_−∞^{∞ 𝑥(𝜏)}_𝑡−𝜏𝑑𝜏 (3.4)

Hence, the amplitude function (absolute values) is obtained when the analytic function, 𝒜(𝑡), is multiplied by its conjugate. Furthermore, to clearly display the major frequency changed with respect to time, the normalized amplitude function was introduced by calculating the percentage of the amplitude over the frequency domain. The visual representation of the normalized amplitude function is equivalent to the wavelet spectrogram, and the results from the three far distant earthquakes are shown in Figure 3.14, Figure 3.15, and Figure 3.16. It can be observed that the major frequency is slightly varied with time, and therefore it is necessary to consider the moving window technique to reflect the time-varying frequency.

doi:10.6342/NTU201901242 41

In the following sections, both the singular spectrum approach and the wavelet approach were conducted to identify and extract the long-period seismic waves. For the singular spectrum approach, the seismic waveforms were first embedded in the lagged covariance matrix and then decomposed by SVD. After that, the singular vectors and the singular values were grouped two by two since the observation in Section 2.1.3 shows that one periodic component usually contributes two coupled singular values with almost equal magnitude and the corresponding singular vectors always have similar periodic property. Generally, the reconstructed time series from the major components contain the long-period seismic waves and can be used to identify the time-varying frequency via FFT and peak-picking technique.

For the wavelet approach, the concept of entropy was introduced to detect the existence of long-period seismic waves. The entropy indicates the amount of information stored in the time series, meaning that the higher entropy represents more information is stored in the time series and vice-versa. The wavelet packet node energy entropy is defined as:

𝐸_𝑖 = (∑ |𝑤_{𝑘 𝑖,𝑘}|^𝑝)

1

𝑝 (3.5)

where 𝑤_𝑖 can be wavelet packet coefficients, 𝑊_𝑗,𝑖 and 𝑆_𝑗,𝑖, or a reconstructed time series, 𝑥̃_𝑖^𝑗(𝑡), from the ith wavelet packet component and 𝑘 stands for the element number of 𝑤_𝑖. The wavelet packet node energy entropy which is actually a p-norm in linear algebra represents the signal energy stored in the particular frequency band and can be used to identify the primary component of the time series in this study. Physically, WPT and Equation (3.5) illustrate that the total signal energy can be decomposed into a set of wavelet packet node energy entropy that corresponds to different frequency bands. The feasibility of applying WPT to the classification of vibration signals has been investigated and the wavelet packet node energy entropy can provide a more robust feature for the classification than using the wavelet packet coefficients directly (Yen and Lin, 2000).

3.3.2 Identification of Long-period Seismic Waves

The same three far distant earthquakes recorded from the B550 Station (from the NCREE network) and the other one far distant earthquake (the Great East Japan earthquake) recorded from the RLNB station (from the BATS network) were selected to identify the time-varying frequency.

Figure 3.6 shows the seismic waveforms collected during the Indian Ocean earthquake. Figure 3.17, Figure 3.18, and Figure 3.19 show the seismic waveforms collected during the other three far distant earthquakes. Again, considering that the P-waves was the fastest of seismic waves, only UD direction was conducted.

doi:10.6342/NTU201901242 42

Singular Spectrum Approach

With a window length of 50 seconds and an appended data length of 20 seconds, moving window SSA was used to extract the primary components of the records starting from the initial P-waves arrival. Consequently, the window length and appended data length were 5000 and 200 points for the B550 station, respectively; the ones were 2500 and 100 points for the RLNB station because of the different sampling rate. The Fourier spectrum and the simple peak picking technique is utilized to identify the time-varying frequency of the four far distant earthquakes, as shown in Figure 3.20. It was observed that the seismic events of the Great East Japan earthquake and the Indian Ocean earthquake show significant long-period seismic waves (larger than 10 seconds) before the S-waves arrival (around 800 seconds). On the contrary, the Lushan (China) earthquake contains no such long-period seismic waves. Moreover, the identified frequency of the Samar (Philippine) earthquake before the S-waves arrival was slightly lower than 0.1 Hz. It was concluded that through moving window SSA, it is possible to detect the content of major frequency of recorded seismic waveforms from the broadband stations, even the major frequency is varied with time. Besides, by applying moving window SSA, only limited computation is needed to extract the major frequency from each time window although embedding and decomposing a large Hankel matrix for the long-period seismic waves may become a time-consumed process.

To examine the major frequency locally measured in the high-tech fabs, the same three far distant earthquakes recorded from the basement of the high-tech fabs were analyzed by moving window SSA, too. These records were originally acceleration collected with a sampling rate of 200 Hz. They were integrated using trapezoidal numerical integration and processed with baseline correction in order to convert to velocity so that the seismic waveforms are comparable with the ones collected from the NCREE and BATS networks. The baseline correction was conducted by Butterworth high-pass filter with a cutoff frequency of 0.01 Hz. Figure 3.21 shows the seismic waveforms from the basement of the high-tech fabs and the B550 station during the Indian Ocean earthquake. Admittedly, only UD direction of the far distant earthquakes was analyzed, as shown in Figure 3.22. Although the integrated velocity was divergent after 1100 seconds for the Lushan (China) earthquake, the result was still informative because the ground motion mainly excited before 1100 seconds. Moving window SSA with the window length of 50 seconds and the appended data length of 20 seconds was used to extract the primary components starting from the initial P-waves arrival. Figure 3.23 displays the time-varying frequency identified by applying FFT and the peak picking technique to the primary components. The observation shows that not only the seismic waveforms but also the identified frequency is consistent between the basement of the high-tech fabs and the station of the NCREE

doi:10.6342/NTU201901242 43

network. This further demonstrates that the recorded waveforms from the different networks are all in a good agreement, as mentioned in Section 3.1.4.

Wavelet Approach

An alternative way to detect the existence of long-period seismic waves from the far distant earthquakes is the wavelet approach (Huang et al., 2015; Huang et al., 2016). The same four records (in UD direction) collected from the NCREE and BATS networks were selected to provide a comparison between the singular spectrum approach and the wavelet approach. To extract the long-period seismic waves from the seismic waveforms, the moving window technique was also adopted.

The window length was 25.6 seconds (2560 points) and data was appended every 5.12 seconds (512 points) for the B550 station; on the other hand, the window length was 25.6 seconds (512 points) and data was appended every 5.6 seconds (112 points) for the RLNB station. Moreover, the wavelet function was Biorthogonal 6.8 and the reference levels of decomposition in WPT were selected as 6 and 4 so that the frequency bands in the first decomposed components ranged approximately from 0 Hz to 0.781 Hz and 0 Hz to 0.625 Hz for the B550 and RLNB station, respectively. The wavelet packet node energy entropy was evaluated using the reconstructed time series with the Euclidean norm where 𝑝 is equal to 2 in Equation (3.5). The wavelet packet node energy entropy from the first decomposed components was then divided by the summation of all the decomposed components to get the percentage of the energy entropy stored in the first components. The percentage of the energy entropy is shown in Figure 3.24. It reveals that the four far distant earthquakes possessed a large portion of long-period components. The results of the other two far distant earthquakes on April 19, 2013, and November 15, 2014, collected from the B550 station were also presented in Figure 3.24 as a comparison. These two recorded seismic waveforms had a much small portion of long-period components as compared to the others and eventually produced no impact on the high-tech fabs because the air mounts filtered out most of the excitation.

Although the contribution of the P-waves to the vertical direction is theoretically largest compared to the other seismic waves for the earthquakes with short distances from the epicenter, because of various transmission paths, it is still questionable if the vertical component of the far distant earthquakes is the fastest one or not. Especially, the distance for the Great East Japan earthquake is 2583 kilometers, the distance for the Indian Ocean earthquake is 3831 kilometers, the distance for the Lushan (China) earthquake is 1876 kilometers, and the distance for the Samar (Philippine) earthquake is still 1621 kilometers. To examine the angle of inclination, the wavelet approach was again applied to every direction of the six far distant earthquakes and the percentage of the energy entropy was evaluated to investigate which direction was the fastest one that reached the

doi:10.6342/NTU201901242 44

percentage to 80%. The results show that only the UD direction of the Great East Japan earthquake and the Lushan (China) earthquake cannot provide the fastest warning. The NS direction of the Great East Japan earthquake and EW direction of the Lushan (China) earthquake may give earlier warning;

however, the difference was trivial, as shown in Figure 3.25.

Brief Summary of Different Approaches

In this section, moving windows SSA was implemented to extract the primary components of the far distant earthquakes. For the identification of the time-varying frequency, the singular spectrum approach takes advantage of the Fourier spectrum and the peak picking technique. Alternatively, moving window WPT was carried out to directly decompose the low frequency seismic waves (approximately below 0.7 Hz). The percentage of the energy entropy is exploited in the wavelet approach to indicate the portion of long-period components. Both the two approaches have very similar results and can be used to detect the long-period seismic waves.

However, the computational complexity of these approaches is completely different. With a window length of 50 seconds and an appended data length of 20 seconds in moving window SSA and a window length of 25.6 seconds and an appended data length of 5.12 seconds in moving window WPT, the computation time was evaluated under the same software and hardware environment, as shown in Figure 3.26. The hardware of the test computer was Intel® Core™ i7-6700K at 4.00 GHz with a memory of 24 GB; the software was Windows 7 (Ultimate) SP1 and Matlab 2010b (7.11.0.584).

The average computation time of the singular spectrum approach and the wavelet approach was 18.6 seconds and 0.07 seconds, respectively. The results clearly illustrated that the wavelet approach took less than 0.1 seconds but the singular spectrum approach spent almost 20 seconds to give a single result. Actually, the only computation time larger than 0.1 seconds in the wavelet approach was the first step because the preparation spends much time than the main tasks (the first step took 0.234 seconds). It is concluded that the wavelet approach is more capable of real-time applications although the singular spectrum approach possesses some adaptability.