Using a Fractal Analysis and Polarization Method for Phase Identification in Three-Component Seismograms

This study presents the automatic P-wave and S-wave arrivals picking algorithm which is essentially based on the fractal dimension and polarized method. With an estimate of the spectral exponent γ in a 1/f process, an interval that indicates the preferred intersection containing both noise and the P-wave is well-detected by corresponding to the minimum absolute spectral exponent γ value along the trace. Based on the different properties of background noise and deterministic signal, the fractal dimension technique can detect the position of the P-wave. The place where the fractal dimension value changes suddenly within the intersection interval indicates the location of arrival of the P-wave. Testing that adds various levels of noise to the seismic signal shows the method can prove able to tolerate noise to a signal-to-noise (S/N) ratio 1.5. Based on the P-wave arrival, the polarized P-wave could be detected by a genetic algorithm (GA) with the strength of polarization and phase difference between the vertical and horizontal components as constraints. Using the first arrival phase as the basis phase, this study combines a polarization filter including rectilinearity functions, linear polarization, phase difference and directionality with GA to detect polarized S-wave of seismograms. Finally, the technique was applied to teleseismic data and near-field motion to verify the accuracy and wide applicability of this method. To conclude, this proposed method, an efficient and brand-new method of associating signal processing technique with physical wave motion properties, may be of importance in finding P-wave and S-wave phase arrivals accurately using three-component seismograms.


INTRODUCTION
Detecting and accurately identifying the first arriving P-wave and S-wave are important in various areas of seismology, including source mechanism analysis, seismic refraction, tomographic studies, event identifications and hypocenter determinations.In a recent study, Stefano et al. (2006) applied an automated picking system to seismograms of the Italian national seismic network and found that the derived arrival times were useful in determining the locations and hypocentral depths of earthquakes.The seismic-velocity structure of a subsurface is also determined from the travel times of the first arrival P-waves recorded by well-distributed stations (Hu et al. 1994;Sato et al. 1996).Some criterion for the onset of signal amplitude and features of the seismic waveform can be used to detect the phase arrivals of a seismic signal.A subjective decision for picking the first arrival P-wave is made based on the change in nature of the trace in terms of amplitude and/or frequency and/or phase both within the trace itself and also relative to its neighbors (Boschetti et al. 1996).Zhang et al. (2003) also noted that a P-wave arrival is characterized by a rapid change in amplitude and/or the arrival of high-frequency or broad-band energy.However, manually picking phase arrivals of seismic events does not provide consistent results (Lingxiu and Moon 2000).Therefore, in past two decades, numerous methods have been developed for identifying the seismic phases.For example, the standard short-term average over long-term average ratio (STA/LTA) algorithm was designed as a criterion for detecting arrivals (Allen 1982).Bear and Kradolfer (1987) modified the envelope function to yield a characteristic function that generates an adjustable threshold to help determine the P-wave arrival.Ruud and Husebye (1992) conceived the idea of combining STA/ LTA algorithm with the polarization approach to identify P-waves.Anant and Dowla (1997) employed the wavelet transform and polarization approach to locate the P and S phase arrivals in short-period seismograms, and compared these results with those obtained from the STA/LTA method to show the accuracy of their method.In addition, the autoregressive (AR) technique is also applied to identify the P and S phases.Sleeman and Eck (1999) adopted the AR model and the Akaike information criterion (AIC) to automatically identify P-phases for a rapid warning system.Zhang et al. (2003) associated AIC technique with wavelet transform to solve the limited window problem and identify the P and S arrivals of a local earthquake or regional events.However, the automatic methods do not perform well when the S/N ratio is low and the arrival is not evident.Hence, picking late arriving seismic signals is quite difficult as the unwanted signals arrive prior to the ones being detected (Lingxiu and Moon 2000).The differences of fractal dimensions between seismic signal and background noise make it useful for methods employing the fractal dimension to distinguish between seismic signal and noise.Bochetti et al. (1996) calculated and analyzed the fractal dimension variations of the seismic trace to identify the seismic first arrival, achieving good agreements with other methods.Tosi et al. (1999) exploited the fact that the fractal dimension of a seismic signal differs from that of random background noise, and well identify the seismic signal from the testing signals which various levels of noises are added to.Investigating a variance fractal dimension (VFD) technique, Lingxiu and Moon (2000) applied their method to real data sets to detect seismic refraction signals from background noise.After analyzing the fractal dimensions of seismic amplitude, phase and instantaneous frequency, Nath and Dewangan (2002) concluded that the fractal dimension analysis can be applied to stacked seismic data for estimating the positions of reflectors.Nevertheless, the S-wave or surface wave arrival of a deterministic part of seismic signal can not be detected directly using only the fractal dimension technique.Based on the discussion above, integrating and taking the advantage of different techniques to develop a composite method for detecting the S-wave phase automatically is a worthy endeavor.
The study proposes a new method that can be used to identify the first-arriving P-and S-wave phase.The proposed method contains two parts.One part of the method combines the 1/f spectral exponent parameter estimation with the fractal dimension to locate the first arrival wave.This method successfully overcomes the two weaknesses of the approach proposed by Boschetti et al. (1996).The first weakness is the manual selection of the working region of the trace that contains the first-arriving wave.When large amounts of data are to be processed, the automatic detection of proper interval is important.The second weakness of Boschetti's approach is that the location of the first Pwave arrival in the seismic trace must almost be known in advance.Therefore, an appropriate working region could be selected properly before applying the fractal algorithm to the selected region of the trace for picking first P-wave arrival.The method developed in this study requires no prior knowledge.
The second part of the proposed method is to modify the polarization parameters (Shieh 1996) with the directionality function (Park et al. 2004) and to associate them with a genetic algorithm (GA) to detect the S-wave arrival.The proposed method can detect the P-wave and S-wave arrivals directly from three components recorded at one station and does not need a reference wave used in the polarized correlation method (Shieh 1996).Compared to the results in Park et al. (2004), the proposed method can match the estimated arrival times of the phase for teleseismic data relatively well and improve the selections of P-wave and S-wave arrivals.

METHODOLOGY
The 1/f family of statically self-similar random processes x(t) measure power spectra that obey a power law of the form where ( ) sx ~ is the power spectra of x(t) over frequency ~, v is a finite nonzero constant, and γ is the so-called spectral exponent.The wavelet transform of a process x(t) is where xn m are the wavelet transform coefficients of the signal x(t), and n and m are the translation and dilation indices, respectively.( ) } represents the normalized dyadic dilation and integer translations of the mother wavelet ( ) t } , Taking the variance of Eq. ( 2) and the base-2 logarithm of both sides yields a straight line whose slope is the estimated γ parameter (Wornell 1993) where the 2 v is a positive real constant that is proportional to x 2 v .If the process x(t) is corrupted by additive white noise g(t), the process is represented as and the parameter γ can be estimated by Baykut et al. (2007): where r m v denotes the variance of the wavelet coefficients rn m of process r(t).The variables γ and 2 v can be derived using the base-2 logarithm of both sides of Eq. ( 6) and fitting the new linear equation using the least-square method.Based on the relationship between spectral exponent γ and fractal dimension D (Wornell 1996) the spectral exponent γ can be used as a quantitative measure of the window that moves along, and includes various parts of, the seismic trace.In addition to the above, Padhy (2004) suggested that background noise has a higher fractal dimension than the seismic signals.Accordingly, from Eq. ( 7), the γ value of a moving window that contains only background noise is lower than that of a window that contains only the seismic signals.Choosing the proper moving window, the minimum of the γ values estimated by Eq. ( 6) will correspond to the interval that contains both background noise and the first arrival.To reduce computation time and ensure that the window encloses both the noise and seismic signal, this selects a moving window length of approximate 1/8 of the seismic signal sampling points.The moving step size of the window across the trace used in this study is specified as 1/5 of the length of the moving window.
After the interval with minimum γ value is found, the first arrival can be identified by applying the fractal dimension technique to the interval.This study uses the rescaled range analytical method of Hurst et al. (1965) to measure the fractal dimension.The Hurst exponent, H, is related to the fractal dimension by The place where the fractal dimension suddenly changes indicates the correct location of first-arriving P-wave (Boschetti et al. 1996).This means that the first-arriving P-wave can be detected by observing the varieties of fractal dimensions of the interval.
Based on the first arrival, the polarized correlation terms offered by Shieh (1996) can be used to locate the polarized P-wave.The following section describes the main steps in this method (Shieh 1996).First, use the Hilbert transform of the three components seismic signal selected at the window range of T seconds, called the P-window in this study, to obtain the analytic triaxial signals.The station coordinate (z, x, y) stands for vertical, east-west and south-north, respectively.As the signals in the P-window are complete polarization, only one eigenvalue of the coherency matrix whose elements are values of the covariance between the analytic triaxial signals of two components of the three-component time series within the P-window exits and the corresponding eigenvector has the form of cos sin cos sin where z is the phase difference between the vertical and horizontal components, az , ax and ay are the normalized direction projections.The phase difference ref z is calculated by Second, based on the eigenvalues of the three-component seismic signals in the selected P-window, write the strength of polarization P ref and rectilinearity R ref as where λ 1 , λ 2 and λ 3 are the eigenvalues of the coherency matrix, in order from maximum to minimum, and e is the ellipticity.Determining the best P-window with range T seconds is an optimization problem whose solution ensures that the signal enclosed in the P-window is polarized.After associating the P-wave arrival with GA, the proper range T seconds of P-window containing the polarized P-wave can be found using the Eq. ( 10) and maximizing the Eq. ( 11) as constraints in GA.If the T seconds is ascertained, the vector, e T i ^h, denoted as the eigenvector associated with the largest eigenvalue λ 1 , can be applied to construct a polarized filter for detecting an S-wave.In other words, the first arrival phase belonging to a P phase is the basis phase for detecting the S phase.
The next step constructs the polarized S-wave filter and applies the GA to determine two significant parameters: the location of the polarized S-wave and the window with a range of T s seconds.This window contains the polarized S-wave, and is called an S-window in this study.Shieh's (1996) method requires that several polarized correlation terms should be defined.The approach in this study, however, slightly changes the contents of the polarized correlation terms to construct the polarized S-wave filter.We redefine the several time-dependent polarized terms in the S-window as follows where t z^h, P t ^h and R t ^h are phase difference, rectilinearity and strength of polarization.The parameters, t z^h, P(t) and R(t), are similar to Eqs. ( 10), (11), and ( 12), but they are calculated in the S-window.Equation ( 16), C t c ^h, can suppress a circular-polarized wave when the correlating P-wave and circular polarization results in φ c = 90°; consequently, C t 0 c = ^h .If the correlated waves are P-waves and other nonlinear waves, the value of C t c ^h will be smaller than 1 (Shieh 1996).e T b s ^h is the eigenvector associated with the largest eigenvalue of the coherency matrix which is constituted by means of the covariance between the analytic triaxial signal of two components of three-component time series within the S-window.Equation ( 17) represents the inner product of two eigenvectors, where 7 and ) stand for the inner product and complex conjugate, respectively (Park et al. 2004).Finally, we define the polarized S-wave filter as When the signals contained in the S-window are completely polarized with the same angle of incidence as the polarized P-wave applied to calculate the time dependent polarized terms, φ c = 180, R c = 1, P c = 1; therefore, C c = 1 and I c = 1, due to the particle motions between polarized P-wave and S-wave are perpendicular.In other words, the position with the maximum value of T c is the desired solution which indicates the S-wave arrival.Determining the S-window range T s seconds, and the proper location of the S-window are also an optimal problem.This study employs GA with maximum S-filter as a constraint to solve the problems.The GA is an adaptive searching procedure as well as a powerful tool to solve problems associated with many unknown parameters (Jin et al. 2000).There have been some previous applications of GA used in source process and geophysical pa-rameters inversions, for example, Šílený (1998), Ramillien (2001), Aguirre and Irikura (2003), and Liao and Huang (2008).Based on the calculations above, the proposed technique can clearly identify the P-wave and S-wave arrivals and overcome the uncertainties of determining the window range of the S-window in Park et al. (2004).

Simulation Results
To confirm the validities of the method, Fig. 1a superposes a random signal which is constructed as a synthesized 1/f 2 process with white noises having a signal to noise ratio (S/N ratio) of 5. Figure 1b displays the power spectral density (PSD) of the simulated signal.Clearly, noise affects the PSD such that decay with a frequency of 5 ~ 100 Hz does not follow the 1/f 2 rule; rather, it is close to a constant.Figure 1c plots the base-2 logarithm of the variances of the wavelet coefficients versus the scales m in Eq. ( 4).The plot of the compound signal as a broken line around m, which is equal to 4, reveals that the lower scales are affected more than higher scales.Figure 1c (solid lines) demonstrates two slopes of the curve due to white noise corruption.The slope of the line (dashed) that best fits the curve in Fig. 1c is approximately 1.41, which is not an accurate γ value.When the formula ( 5) is applied to the signal of Fig. 1a, the correct estimation of slope is approximately 2.01, as Fig. 1d indicates.
This performance shows that the method (Baykut et al. 2007) provides very good estimates of the spectral exponent γ value of the signal with heavy noise, suggesting that the method is useful for determining the optimal γ value of the signal to investigate the hidden characteristics in the signal further.

Picking the P-Wave Arrival
This section investigates the accuracy of the proposed method in picking the P-wave arrival, and discusses the effects affected by different levels of noise.This test analyzes a teleseismic wave recorded at station INCN from the earthquake (M w = 5.9) that occurred on 22 October 1999 in Chia-Yi, Taiwan and displays the results in Fig. 2. A window is then moved across the seismic trace and the γ value of a different part of the seismic trace within the window is calculated.Figure 2d displays the γ values along the seismic trace.Figure 2a depicts the region that corresponds to the minimum γ value as a box bordered with dashed lines.To clearly reveal the signal enclosed in the region, Fig. 2b enlarges the diagram of the region.This figure shows that the region contains both the noise and the deterministic part of seismic signal.Therefore, in view of the difference between deterministic signal and   background noise in terms of fractal dimension, Fig. 2c calculates the fractal dimensions of the signal displayed in Fig. 2b to determine the P-wave arrival, and displays the results.Furthermore, the fractal dimension curve in Fig. 2c shows that the fractal dimensions of noise are higher than those of a deterministic seismic signal.At a critical point, the fractal dimension changes rapidly with the seismic attributes in the seismogram.This occurs close to the transition region from noise to the first arrival.Figure 2c plots the variations in the fractal dimension, indicating the sudden change with a square.The P-wave arrival can be identified by the location of this square (triangles in Fig. 2b).The test described above produces two important findings.The first one is that the proposed method can identify automatically the correct interval that contains both background noise and seismic signal.The second is that the P-wave arrival can be found precisely by analyzing the variation in the fractal dimensions of the signal within the interval related to the minimum γ value.
To evaluate the robustness and the effectiveness of the proposed method, this study again uses the same teleseismic wave in the test above as an example, but adds an increased amount of noise to the wave signal.The purpose of this test is to estimate the maximum amount of noise the proposed method could tolerate.This is quite important in the presence of large amounts of noise because the real seismic data sets with varying S/N ratios recorded at worldwide stations are foreseeable when an earthquake occurs.However, in a situation of very low S/N ratio, the detection of first-arriving P-wave may be affected seriously.Figure 3 plots these re-  sults.Figure 3a represents the seismic signals with different level noises: S/N ratios of 20, 10, 5, 3 and 1.5.The region of each seismic signal that corresponds to the minimum γ value is plotted as a box bordered with dashed lines in Fig. 3a, respectively.The positions of the intervals with the smallest γ values estimated in different noise circumstances are not consistent, but all of them contain the first arrival and background noise well.The γ values along the seismic traces (Fig. 3a) are demonstrated in Fig. 3b. Figure 3c is the fractal dimensions of the signal displayed in Fig. 3d to determine the P-wave arrival.The location of a square box in Fig. 3c symbolizes the sudden change of the fractal dimensions, indicating the P-wave arrival.Figure 3d shows the enlarged diagrams of the bounded regions shown in Fig. 3a, in which red triangles indicate the P-wave arrival.This test shows that the proposed method can locate the preferred interval and identify the P-wave arrival well, even though various amounts of noise were added.However, the background noise may reduce the accuracy of the identification of the P-wave when S/N is less than 1.5.

Application 1
The study employs the proposed method to analyze the earthquake (M L = 6.4) recordings obtained on 22 October 1999 at Chia-Yi, Taiwan.This experiment uses data from ten teleseismic stations (Fig. 4) uniformly distributed around the epicenter of the Chia-Yi earthquake.Figures 5  and 6 depict the results.In these figures, a rectangle (dashed line) and a triangle indicate region detected by the 1/f process and the P-wave arrival, respectively.The results in Figs. 5 and 6 clearly show that the method, when applied to the teleseismic wave, correctly and automatically identifies the P-wave arrival.It is also effective when applied to noisy data, such as that obtained from stations YSS, PMG and TLY, easily locating the first arrival and regardless of the particular characteristics of the signal such as frequency or amplitude.Based on the effectiveness of determining the P arrival, the S arrival can be well detected by combining S-filter with GA.

Application 2
This section applies the proposed method to teleseismic waves, noisy seismic data and near-source station recording to identify the P-wave and S-wave arrivals, and further demonstrates the excellent performance of the proposed method.To provide an accurate comparison with the theoretical arrival times calculated by Park et al. (2004), the data set consists of the teleseismic event of the 21 September 1999 Chi-Chi earthquake with epicentral distance range of 30.89° (station TLY) and 50.07° (station CTAO).Figures 7 and 8 demonstrate the results of total polarized correlations, R c (t), P c (t), C c (t), I c (t) and T c (t) from top to bottom, of the teleseismic waves recorded at stations TLY and CTAO.The S-filter, T c (t), is normalized up to 1.At the maximum of T c (t), the locations of S-wave arrivals can be located clearly about 338 s for the TLY station and 470 s for the CTAO station.These results agree quite well with the results in Park et al. (2004), which referred to the Jeffreys-Bullen Travel Time Tables to identify the seismic wave phases.These results also show that a shorter length of Swindow produces denser oscillations in polarized correlations appearing in Fig. 8. Conversely, a longer S-window length has a smoothing effect on the polarized correlations in Fig. 7.The proposed method accurately picks the S-wave arrival in Fig. 7 without other high value in T c (t) that might interfere with picking the first S arrival.Furthermore, note that the high values from 510 to 590 s of T c (t) in Fig. 8 indicate the surface waves contain significant S-wave energy after the first S arrival.Figures 9 and 10 display the three-component seismograms recorded at stations TLY and CTAO.The main achievement of this study is using the P-wave in the P-window as a reference wave to carry out the polarized correlation analysis, and subsequently, identifying the polarized S-wave arrival.Figures 9 and 10 show two square windows calculated by the proposed method.These squares represent the P-window and S-window, and reveal the P and S first arrivals respectively.Figure 9 clearly shows that it would be difficult to locate the arrival times of the phases using only the variations of frequency contents or waveform amplitudes.The method proposed in this study can directly and accurately identify the first-arriving P and S phases without using other seismic data as a reference used in Shieh's method (1996).Figure 10 shows that first-    arriving P and S phases are identified clearly appearing at proper time, as the results in Park et al. (2004).The results of Park et al. (2004) are carried out using the directionality and rectilinearity functions; however, the accuracy for first S arrival picking is contaminated by other S-phases, so that the method is confined to the quality of recorded data and requires the Jeffreys-Bullen Travel Time Tables to identify the seismic phases.In this study, the results of Figs. 9 and 10 not only agree with those of Park et al. (2004)(dashed lines in Figs.9 and 10), but also illustrate that the proposed method refines the method offered by Shieh (1996) and overcomes the limitation of Park et al. (2004).
In order to demonstrate the validity of the proposed method for noisy seismic data, we apply the proposed method here to analyze the observed seismic data of the Yilan earthquake (M L = 5.1) in Taiwan on 23 January 1996.The seismic data used in this study is recorded at CHY station with the epicentral distance of 206.2 km.In the analysis of the seismic data, the arrival times of P and S waves are determined and compared with the estimated arrival time offered in the database of CWB (Central Weather Bureau). Figure 11 displays the correlation coefficients and S-filter.In Fig. 11, it is apparent that the arrival time of S-wave corresponding to the maximum of T c (t) is clearly located at 57.36 s, without other higher value in T c (t) that may in-terfere with the detection.Figure 12 depicts the result of phase identification analysis for the seismic data.As shown in this figure, the seismograms contain significant low-frequency noises which may interfere with the seismic phase identification.However, from the comparisons between the results derived from the proposed method (solid line) and the estimated arrival times (dashed line) of CWB in Fig. 12, we can find that our results are in complete agreement with the estimations of CWB, and the differences between the two estimations are only within 0.3 s.Accordingly, the findings suggest that the proposed method can significantly contribute to effectiveness of picking the Pand S-wave arrival times of the noisy seismic data.
Finally, this study applies the method to analyze the near-field strong motion.This experiment uses data from the earthquake on 22 October 1999 in Chia-Yi, Taiwan as recorded by testing Station TCU065, located in the central area of Taiwan.Figure 13 shows the correlation coefficients and S-filter, revealing that the maximum of T c (t) at about 27.2 s, which is the location of first S-wave arrival.However, the variation of T c (t) in Fig. 13 is much higher than those in teleseismic recordings (Figs. 7 and 8).This may be because the seismic phases recorded at near-field are more complicated than at far-field.In the three portions of T c (t) in Fig. 13, the first portion that is before 27 s consists of   low value of T c (t), indicating that the waves between first P arrival and 27 s contain fewer S-wave components.In the second portion, between 27 s and S-wave arrival, the high values of T c (t) are contributed by the converted S waves (P to S).The last portion, after the first S arrival, may be composed of surface waves, and the high values of T c (t) in this section indicate that the surface waves include quantities of S-wave energy.In addition, the value of C c (t) in Fig. 13 corresponding to S-wave arrival time is not equal to 1.This may be because the angle of incidence of S-wave is affected by anisotropic media to make the value of φ c in Eq. ( 13) unequal to 180; consequently C c (t) < 1.Nonetheless, the effect of C c (t) does not affect the picking location of S-wave arrival in T c (t). Figure14 shows the three-component seismograms collected by station TCU065.Two square windows, the P-window and S-window calculated by the proposed method, indicate the first P-wave and S-wave arrivals, respectively.As Zhang et al. (2003) noted, a P-wave arrival is characterized by a rapid change in amplitude or the arrival of high-frequency or broad-band energy.Viewed in this light, the first P-wave arrival located at about 20.5 s in Fig. 14, which changes both in amplitude and frequency, provides a reasonable result.However, the S phase arrival is significantly more difficult to locate than the P phase arrival (Anant and Dowla 1997).The S-wave arrival in Fig. 14 coincides with the empirical rule that the first S arrival locates at the place where a rapid change in wave amplitude or in the arrival of lower frequency.Therefore, it is likely that the results are correct and the method is effective in identifying the P and S arrivals.

CONCLUSION
This study employs the 1/f process, fractal dimension, polarization method, and GA to obtain the appropriate windows that captures both the first P and S arrivals.One part of the method developed in this study is based on the wavelet coefficients of the signals that enclosed in a window moving across the trace.Taking the variance of the wavelet coefficients and base-2 logarithm of the variance, we can calculate the spectral exponent γ parameter.The interval that corresponds to the minimum value of the spectral exponents, containing both background noise and seismic signal, can be applied to detect the first P arrival.Theoretical testing reveals that the method can confirm the preferred interval from noisy data, and accurately identify the first-arriving wave when combined with the fractal-based picking technique.This study employs the proposed method to identify the first P arrivals of a set of teleseismic recordings.The results conclusively demonstrate that the first arrivals of ten teleseismic recordings are detected well.On the basis of the first arrival, the polarized P wave is located by combining GA with polarization method.Due to different particle motions of the polarized waves, this study offers the S-filter T c (t) composed of modifying correlations terms of Shieh (1996) to identify the first S-wave arrival.The position with the maximum of T c (t) indicates the right S-wave arrival.Two teleseismic recordings used in Park et al. (2004) are applied to prove the accuracy of the proposed method.The results derived by the proposed method are entirely consistent with those in previous study without reference to the Jeffreys-Bullen Travel Time Tables.In addition, we also analyze the strong motion recordings and noisy seismic data by employing the proposed method and get reasonable results in picking fist-arriving P and S waves.The results of this study clearly support the effectiveness of the proposed method.Accordingly, the technique in this study improves the method by Shieh (1996), which requires the three-component seismograms of one earthquake as reference waves to detect the P and S arrivals of another earthquake.This method also solves the imperfections presented by Park et al. (2004), in that the problem occurs as the first S arrival is interfered by another S phase, and the time window length of covariance matrix is determined unsteadily.Therefore, the scheme is a valuable tool for determining the first P and S arrivals and has further applications in seismology.

Fig. 1 .
Fig. 1.The theoretical test of estimating spectral exponent parameter of the 1/f process with white background noise.(a) is the simulated signal obeying the 1/f 2 process with added white noise with a S/N ratio of 5; (b) represents the power spectral density (PSD) of the simulated signal; (c) is the base-2 logarithm of the variances of the wavelet coefficients versus the scales m, where the dashed line represents the best linear fitting between scales and logarithm values; (d) shows the best linear fitting between scales m and logarithm values derived from the modified method byBaykut et al.(2007).This figure shows that the noise effect has been eliminated, and provides the correct calculation of γ value.

Fig. 2 .
Fig. 2. Test applying the proposed method to a teleseismic wave without adding white noise.(a) waveform and the dashed-line box in the waveform represents the optimal interval corresponding the minimum of γ value.Triangles indicate the first P-wave arrival position; (b) signal enclosed in the dashed-line box.Triangles denote the position of the first-arriving wave; (c) fractal dimension curve of the signal enclosed in the dashed-line box.The square indicates the sudden variation of the fractal dimension curve, and the position of this square represents the first P-wave arrival; (d) spectral exponents of the moving window across the wave recording.

Fig. 3 .
Fig. 3. Test applying the proposed method to a teleseismic wave while adding white noise with S/N ratios of 20, 10, 5, 3 and 1.5.(a) Seismograms with different levels of noise were added.The dashed-line boxes in the waveforms are suitable intervals corresponding to the minimum γ values under different noise conditions.Triangles indicate the first P-wave arrival position; (b) spectral exponents, γ values, of the moving window across the wave recording; (c) fractal dimension curve of the signal enclosed in the dash-line box.The point where the fractal dimension curve changes is indicated by a square box; and (d) enlarged diagrams of the dashed-line boxes of waveforms in Fig. 3a.Triangles indicate the first P-wave arrival position.

Fig. 4 .
Fig. 4. Location of the epicenter of the 22 October 1999 Chia-Yi, Taiwan, earthquake (star) and the ten teleseismic stations (triangles) used in this study.

Fig. 6 .
Fig. 6.The results of five seismograms recorded at the stations YSS, MAJO, PMG, CTAO and CHTO.The dashed-line box is detected by the 1/f process, and the triangle represents the first P arrival.

Fig. 7 .
Fig. 7.The total polarized correlations, R c (t), P c (t), C c (t), I c (t) and T c (t) of the teleseismic waves recorded at stations TLY.The position with the maximum of T c (t) indicates the S-wave arrival.

Fig. 9 .
Fig. 9. Three-component waveforms collected by station TLY.The first window contains the first P arrival, and the second contains the first S arrival.The dashed line indicates the arrival times of P and S waves estimated by Park et al. (2004).

Fig. 10 .
Fig. 10.Three-component waveforms collected by station CTAO.The first window contains the first P arrival, and the second contains the first S arrival.The dashed line indicates the arrival times of P and S waves estimated by Park et al. (2004).

Fig. 11 .
Fig. 11.The total polarized correlations, R c (t), P c (t), C c (t), I c (t) and T c (t) of the noisy waves recorded at station CHY.The position with the maximum T c (t) indicates the S-wave arrival.

Fig. 12 .
Fig. 12. Three-component waveforms collected by station CHY.The first window contains the first P arrival, and the second contains the first S arrival.The dashed line indicates the estimated arrival times of P and S waves by CWB.

Fig. 13 .
Fig. 13.The total polarized correlations, R c (t), P c (t), C c (t), I c (t) and T c (t) of the strong motion waves recorded at stations TCU065.The position with the maximum T c (t) indicates the S-wave arrival.The T c (t) in this figure is more complicated than that of teleseismic recordings displayed in Figs.7 and 8.

Fig. 14 .
Fig. 14.Three-component strong motion waveforms collected by station TCU065.The first window contains the first P arrival and the second contains the first S arrival.