Intrinsic and attenuative dispersion characteristics of direct P-waves in and near the source area of the 1999 M_W7.6 Chi-Chi, Taiwan, earthquake before and after the mainshock

1.   Introduction

The medium of the Earth has inhomogeneous, anelastic and anisotropic characteristics. Dispersion is known to fall into two basic classes: intrinsic, which is based on anelasticity; and scattering, which is based on local wavelength-scale variations in the rock formation (Parra et al, 1999). Intrinsic dispersion is a local property of the rock. Scattering dispersion is a property of a neighborhood of rocks, and includes the effects of re- flections, refractions, and the law requiring continuity of displacement. The basic parameter that is used to describe the attenuation of seismic waves is the dimensionless quantity named quality factor Q. There are several reasons that detailed knowledge of attenuation mechanism of the medium is required. First, measurements of Q constitute a fundamental tool in the evaluation of the information on the component and property of the rocks, as well as temperature distribution. In volcanic areas, low Q values are usually associated with fractured systems fluid-filled (Zucca et al, 1994) or to magmatic bodies (Sanders and Nixon, 1995). Second, measurements of Q are important as indicators of macro-scale Earth heterogeneities, not easily amenable to study using only seismic velocities (Winkler and Nur, 1982; Romanowicz, 1990, 1995; William and Eugene, 2000). Finally, a detailed knowledge of the path effects due to attenuation is required when correcting displacement spectra to estimate source parameters of tectonic earthquakes (Abercrombie, 1997).

For estimating Q from P- and S-waves, many techniques have been developed (Giampiccolo et al, 2003). In summary, there are two kinds of main approaches for studying the wave attenuation: amplitude decay (Liu et al, 1994; Nava et al, 1999; Castro et al, 2003; Horasan and Guney, 2004) and dispersion. Body-wave dispersion has been applied to several studies ranging in scale from direct P-waves, shallow-crustal refraction to global dimensions (Trong and Granet, 1980; Correig and Mitchell, 1989; Cong et al, 2000; Giampiccolo et al, 2002; Liu et al, 2005).

Correig(1991a, b) put forward a method to measure body-wave dispersion in terms of the arrival of narrow band-pass filtered signal. Cong et al (2000) developed this method for measuring attenuative dispersion of direct P-waves. In comparison with the amplitude decay measurement, the measurement of attenuative dispersion analysis has several potential advantages. Measurement of attenuative dispersion can not be severely affected by focusing and defocusing in the regions with laterally complex structure, requires only a single-station recordings of very small events at short distances, what is more, does not require any knowledge of the earthquake focus, can avoid dealing with the source complexity of large earthquakes by using seismic records from small events, and is no problem with signal-generated noise that may contaminate later arrivals. Therefore, if the recorded events cover a wide epicentral distance and depth range it may be possible to map spatial variation of anelasticity.

In the present study, we develop a procedure of Correig(1991a, 1991b) for measuring the dispersion of direct P-waves and applied it to measuring the dispersion of direct P-waves recorded by SML and ALS stations located in and near the source area of the Chi-Chi Taiwan earthquake. The goal of this work is to determine the attenuation characteristics of P-waves in Chelungpu thrust fault region before and after Chi-Chi strong earthquake of magnitude M_W7.6.

2.   Theoretical background for determining quality factor Q(ω) by wave dispersion analysis

Propagation of seismic waves in an attenuating medium can be modeled as the propagation of seismic waves in an elastic medium convolved with an attenuation operator. The attenuation operator can be described by means of a continuous relaxation model (Liu et al, 1976). The attenuation factor of medium is expressed as

where ω is angular frequency, τ₁ and τ₂ are the long and short relaxation times, respectively, τ₁^-1 and τ₂^-1 are the two frequencies at the half amplitude points of the Q^-1(ω) spectrum, and Q_m^-1 is a constant defined by the flat part of the Q^-1(ω) spectrum bounded by τ₁^-1 < ω < τ₂^-1. In order to illustrate the physical meaning of continuous relaxation model parameters, Figure 1 shows simulated Q^-1(ω) based on equation (1).

Figure  1.  The variation of internal friction Q^-1 with frequency for one model defined by τ₁=500 s, τ₂=0.02s and Q_m=50 (open circles denote the frequencies at the half-amplitude points of the spectrum).

DownLoad: Full-Size Img PowerPoint

If the medium is modeled as a linear viscoelastic solid, the phase velocity can be expressed as (Correig, 1991a)

where C_∞ and C(ω) are phase velocity at infinite frequency (perfectly elastic case) and at frequency ω, respectively. Based on the equations (1) and (2), and the following relation between group velocity U(ω) and phase velocity C(ω),

we obtained the analytic expression equation (4)

C_∞ and U(ω) in equation (4) can be rewritten as a function of travel time by substituting C_∞=s/t and U(ω)=s/t(ω) where s, t, t(ω) are travel distance, travel time at infinite frequency and at frequency ω, respectively. If some frequency components of the signal are available, the dispersion can be measured as the difference in arrival of a frequency ω with respect to reference frequency ω_r,

Therefore, the transfer function that relates intrinsic dispersion to attenuation can be written as

Using the observed group velocity delays and travel time of seismic wave, we can invert this equation (6) to obtain Q_m, τ₁ and τ₂ parameters.

3.   Data

Taiwan is located in a tectonically active region. The Philippine Sea Plate is moving northwestward at a rate of approximately 7-8 cm/a (Yu et al, 1997) relative to the Eurasian Plate, creating the Taiwan collision zone. This makes Taiwan Island one of the most active earthquake regions in the world. The M_W7.6 Chi-Chi earthquake occurred in 20 September 1999 (GMT time) with its epicenter at 23.853°N; 120.816°E near Chi-Chi town in Nantou County of western Taiwan. The earthquake produced an approximately 100-km-long surface rupture, mostly along the Chelungpu fault, a low-angle reverse fault in the nearly north-south direction and dipping to the east (Johnson et al, 2001).

In our study we use the data from SML and ALS stations recordings with velocity sensors as a part of the Taiwan Center of the Central Weather Bureau Station Network. SML station is located about 9.2 km ENE of the epicenter of the Chi-Chi earthquake. ALS station is located about 37.9 km south of the epicenter of the Chi-Chi earthquake. SML and ALS stations are equipped with three-component geophones, with a natural period of 1 s and sample signals at a rate of 50 Hz or 100 Hz. The recorder has a flat velocity spectrum between 1 Hz and 100 Hz. For this study we selected 357 and 627 vertical seismograms of direct P-waves with good signal to noise ratios recorded by the SML and ALS stations, respectively. The epicentral distances of the selected earthquakes are mainly less than 30 km, focus depth are mainly less than 30 km, and the duration magnitudes are mainly between 1.5 and 3.0. Figure 2 shows the location of local events during years between 1995 and 2003 and stations used in our study.

Figure  2.  Epicenters of earthquakes (solid circle denotes epicenter of Chi-Chi earthquake and open circles denote earthquake recorded by ALS and SML stations used in this study) and station distribution (solid triangles) (CCWBSN) during 1995-2003.

DownLoad: Full-Size Img PowerPoint

4.   Methods

4.1   Data resampling method

In order to be more easily made the estimation of dispersion, we adopted Newton polynomial insertion value and Gauss filter methods to obtain a sampling signal at 1000 points/s instead of the original sampling points/s. Supposing the original signal sampling rate is f_s points/s, and the resampled signal sampling rate is Kf_s points/s, the Nyquist frequency will be raised to Kf_s/2. In order to remove the noise in the frequency range from f_s/2 to Kf_s/2, we have to perform digitally filter with Gauss window to the data with a prohibitive band of frequency f_s/2. This process can improve the smoothness of the group delay plot (Cong et al, 2000).

4.2   Multiple-filter method of data

We proposed Morlet wavelet filter method to estimate group velocity delay as a function of frequency for the spectral components that make up the direct P-wave. The wavelet band-pass filters based on Morlet wavelet (Morlet et al, 1982; Grossman and Morlet, 1984) have three advantages superior to the usual FIR digital filter: filtering a signal into narrow frequency bands with good frequency response but causing no phase shift, exhibiting no Gibbs phenomenon and suppressing some singularities.

The continuous wavelet transform CWT_x(a, b) of a time signal x(t) can be defined as the sum over all times of the signal multiplied by the scaled and shifted version of the original (mother) wavelet ψ(t):

where a is a real constant and great than 0 (the scale), b is a real parameter (the position). The Morlet mother wavelet consists of a plane wave e^-iω₀t modulated by Gaussian (Kumar, 1995).

We discretize the parameter t into a grid in equation (7). The parameter t is then obtained from their corresponding grid coordinate i, using simple expression of the form t_i=t₀+iΔt (i=1, 2, …, N), where t0 is the lower limit and Δt is the grid unit for the t parameter. In the same way, we may express a and b as a_j=a₀+jΔa (j=1, 2, …, K) and b_l=b₀+lΔb, respectively. The discrete Morlet wavelet transform can be written as

When subscript l of variable b_l ranged from 2M to N+2M, DWT_x(a_j, b_l) was taken for Morlet wavelet transform results.

The scale a and frequency f performed to Morlet wavelet Fourier transform exist following linear relation (Laura and David, 1995):

where ω₀ is non-dimensional frequency in our case ω₀=6, c is a constant.

Therefore, the scale factor has clear physical meaning turned scale variable into frequency variable. DWT_x(a_j, b_l) contains real section Re[DWT_x(a_j, b_l, k)] and imaginary section Im[DWT_x(a_j, b_l, k)] due to Morlet wavelet is complex wavelet. We may calculate the instantaneous phase by means of real and imaginary section. Based on square of the absolution of complex number DWT_x(a_j, b_l), we may get the wavelet energy spectrum value E_x(f_j) corresponding to difference frequency f_j, which may be written as

where

Based on the studied results of Grossman and Morlet (1984) and Mallat (1989), supposing the arbitrary frequency band-pass ranges from j=m to j=m+js which m and js are all constant specified by actual question, then the filtered signal X(t_i) can be retrieved with following equation

The constant r from equation (11) is defined as:

Based on equation (11), we may obtain a series of filter signal with difference band-pass frequency and estimate group velocity delay by means of the arrival time of maxima magnitude from band-pass filtering signals. The selection of parameters for the present study is based on the following criteria: ① the sampling internal Δt is equal to the reciprocal of data sampling points per second. Δa is the same as Δt; ②sampling range of Morlet wavelet is from -M to M which is equal to three times of K in order to embody the full change of wavelet function and avoid the waste of calculation time and information losing; ③ the minimum start scale a₀ is equal to two times of Δt and maximum scale a_k is less than half of time series x(t_i) length. ④ band-pass frequency widths ranged from 0.75 f_j to 1.25 f_j (f_j is center frequency).

4.3   Inversion method of t, Q_m, τ₁ and τ₂ parameters

For all events the arrival times of P and S phases are picked so propagation times t of direct P-waves and epicentral distances can be obtained by means of the Taiwan western speed crust model (Shin and Ho, 1994), arrival time difference between S- and P-waves phase and focus depth results from CCWBSN in Taiwan.

We apply the genetic algorithm (GA) to invert parameters Q_m, τ₁ and τ₂ on the basis of observed group velocity delay and travel time t of direct P-wave from one earthquake event. Detailed discussions on the mechanisms of GA can be found in references (Sambridge and Gallagher, 1993; Chai et al, 1996; Kim et al, 1999; Mackenziea et al, 2001; Aytug and Saydam, 2002; Li et al, 2004; Pezeshk and Zarrabi, 2005). GA has been proved to be a versatile and effective approach for solving optimization problems such as combinatorial and discrete optimization problems, mixed-discrete non-linear optimization problems.

Studied results indicated that various control parameters (such as crossover rate, mutation rate, population size, and the limits for the number of generations) play a crucial role for a successful implementation of a GA. Performance of a GA is significantly affected by these parameter values. However, it is not easy to adapt a search algorithm to a given problem, i.e., to find parameter values that are most appropriate for a given problem. We test several values for each of these parameters in this paper. In the following, we describe how these factors are determined or selected in our implementation of GA. Detailed discussions on the mechanisms of GA can be found in Sambrige and Gallagher (1993).

To use genetic algorithms (GA) for attenuative dispersion measurement we first discretize the entire parameters into a 3-D grid. The parameters, Q_m, τ₁ and τ₂, are then obtained from their corresponding grid coordinate (i_Q, i_τ1, i_τ2), using simple expression of the form Q_m=Q_min+i_QΔQ_m, where Q_min is the lower limit and ΔQ_m is the grid unit for the Q_m parameter. If the number grid points are all powers of 2, the above three parameters are uniquely defined by the three binary numbers (b_Q, b_τ1, b_τ2), where bQ is the binary representation of i_Q, etc. In this paper, initial parameters Q_m, τ₁ and τ₂ selected ranges are [1, 2000], [0, 10000] and [0, 0.5], respectively. The binary value i_Q, i_τ1 and i_τ2 are 32.

We selected following misfit function:

where O_i(k) is observed data, T_i(k) is theory data, M is the number of inversion parameters and S is input population of parameters whose value is equal to 32.

We select the linear scaling scheme suggested by Sambrige and Gallagher (1993) as fitness function. It can be written as

where the constants a and b are determined by the distribution of misfit function φ_k, b=S^-1(φ_max -φ_avg)-1, a=bφ_max, where φ_max, φ_avg are the maximum and mean of misfit values in the current population, respectively. The scaled values, P_kk=1, 2, …, S), constitute a scaled fitness distribution that satisfies and P_k≥0 for all k.

We select crossover rate P_c=0.9 and mutation rate P_m whose values varies in the linear ramp changes from P_ms=0.1 to P_me=0.001. P_m can be written as

where P_m(n) is a function of bit order n, N is maximum bit order.

In our implementation of the GA, the search procedure is terminated when the computation loop numbers exceeds a predetermined numbers limit that is equal to 300.

To demonstrate how the GA may be used for rapid inverting parameters Q_m, τ₁ and τ₂, particularly in the presence of local minima based on the above designed method and parameters, we consider a synthetic example. Supposing Q_m=182.95, τ₁=900, τ₂=0.01, f_r=4 Hz and discretization the parameters f_i(i=1, …, 20) into 20 points in equation (5), we can obtain the theory change values of relative group velocity delay against frequency (Figure 3a, solid points). When τ₁=6 000 instead of 900 and other parameters immovability, the theory change curve of relative group velocity delay against frequency is shown Figure 3a (solid line). The tested results indicate that group velocity delay is very insensitive to changes in parameter τ₁. The random noise value (Figure 3b) whose maxima value is less than 0.002 is added to the synthetic data (Figure 3c, solid points). On the basis of synthetic data, we obtain four groups parameters Q_m, τ₁ and τ₂ by means of following four genetic algoritms, namely, linear ramp from equation (13), exponential ramp from Sambrige and Gallagher (1993), linear fitness function from equation (15) and exponential fitness function from Sambrige and Gallagher (1993) genetic algorithms. In this way we get a more reliable impression of how the algorithms are likely to perform for any initial population. Figure 3d shows the optimum solution with minima stand residual difference (Figure 3c, solid line) and the population average curves (Figures 3d and 3e). We notice that in terms of both minimum inversion error the linear ramp and linear fitness function algorithm is the most successful of the four genetic algorithms. We conclude that 1 the inversion process using a GA results in a good agreement between the theoretical and experimental dispersion curves; 2 The GA part of the procedure is fast, stable and accurate, with several advantages compared with traditional methods.

Figure  3.  The synthetic data and fitting curve obtained by GA, and the change curves of the lowest misfit and average misfit of each population against the times of iteration. (a) A plot of the obtained theoretical dispersion curve using the equation; (b) A plot of random noise dispersion curve; (c) The synthetic data and fitting curve obtained by GA; (d) The lowest misfit is plotted against the times of iteration; (e) The average misfit of each population is plotted against the times of iteration.

DownLoad: Full-Size Img PowerPoint

As an example, Figure 4 shows original and resampling seismograms, observed dispersion of spectral components of the first cycle of direct P-wave and fitting curve of Q_m value etc. parameters. Figure 4a shows the seismogram of an M_L1.7 earthquake occurred at 08:14:44 am of 30 November 1995 (GMT) and the arrival onsets time of direct P- and S-waves recorded by ALS station. The travel time and epicentral distance is 2.38 s and 4.5 km obtained by the time difference between direct S- and P-waves arrival time and the Taiwan western speed crust model (Shin and Ho, 1994), respectively. Figure 4b displays the waveform after resampling with the sampling rate of 1000 samples per second. Figure 4c shows the first cycle signal performed data resampling. The open circles in Figure 4d denotes observed dispersion of spectral components of P-waves, and whose diameters were proportional to observed spectral amplitudes. The dashed line in Figure 4d denotes inversion result by means of GA method.

Figure  4.  Original and resampling seismograms, observed dispersion of spectral components of the first cycle of direct P-wave and fitting curve of Q value etc. parameters. (a) Seismogram of an M1.7 earthquake at 08:14:44 of 30 Nov.1995(GMT) recoded by ALS station with a sampling interval of 0.02 samples/s, the arrow denote the arrival time of direct P- and S-waves; (b) Resampled seismogram for event recoded by ALS station with a sampling interval of 0.001 samples/sec; (c)The first full cycle of the Pg wave in the seismogram of (c); (d) Observed dispersion of spectral components of the first full cycle of the P wave in the seismogram of (b). The open circles indicate dispersion from (b) signals. The diameters of the open circles are proportional to observed spectral amplitudes. The dashed line indicates fitting curve.

DownLoad: Full-Size Img PowerPoint

5.   Intrinsic and attenuative dispersion characteristics of direct P-waves

Figure 5 shows the distribution of Q_m obtained from SML station data records versus epicentral distance, earthquake focus depth and earthquake magnitude. In the Figure 5, we select 457 events whose epicenters mainly distribute within the range of 25 km, and whose 5.6%, 10.4%, 5.9%, 52.9% and 25.2% of events occurred in the 4-9 km, 9-13 km, 13-17 km, 17-25 km and more than 25 km depth ranges, respectively. Figure 6 shows the distribution of Q_m obtained from ALS station data records versus epicentral distance, earthquake focus depth and earthquake magnitude. In the Figure 6, we select 627 events whose epicenters mainly distribute within the range of 25 km, and whose 20.0%, 41.8%, 19.8%, 17.2% of events occurred in the 4-9 km, 9-13 km, 13-17 km and 17-25 km depth ranges, respectively. We can see from Figures 5 and 6 that the distribution of Q_m are characterized by increase with distance and depth, and is clearly independent of magnitude.

Figure  5.  Q_m values obtained from SML station data records varying with distance, depth and magnitude.

DownLoad: Full-Size Img PowerPoint

Figure  6.  Q_m values obtained from ALS station data records varying with distance, depth and magnitude.

DownLoad: Full-Size Img PowerPoint

In order to remove the effects of distance and depth on the Q_m values, we divided the set of Q_m values into some groups by means of different depth range. Shin and Ho (1994) reported that the crust block of the Taiwan westward was considered to consist of seven layers whose depth range were 0-1 km, 1-4 km, 4-9 km, 9-13 km, 13-17 km, 17-25 km and 25-30 km, respectively. As an example, Figures 7 and 8 show Q_m values and its residual characteristics varying with epicentral distance and focal depth obtained from SML and ALS stations data records, respectively. We can see from Figures 7 and 8 that the distribution of Q_m versus epicentral distance still reveal an linear increase tendency of Q_m with increasing the distance in specified depth range (Figures 7a and 8a). But the distribution of Q_m versus specified depth range reveal scatter but no obvious systematic change of Q_m with increasing the depth (Figures 7c and 8c). After the linear regression term in Figures 7a and 8a are removed, we obtain the Q_m residual varying with epicentral distance and focus depth (Figures 7b and 7d, Figures 8b and 8d). It can be clearly seen that Q_m residual with increasing the distance have greater scatter but no systematic linear tendency change. Q_m residual with increasing depth in Figures 7d and 8d show the same change rule as that in Figures 7c and 8c, but Q_m residual exhibit less scatter. Therefore, the above Q_m residual from specified depth range which corresponded to one or two layers of the crust block of the Taiwan westward show no obvious correlation with distance, depth and magnitude.

Figure  7.  Q_m values obtained from SML station data records at depth range 13-25 km, and its residuals characteristics varying with epicenter distance and depth. (a) Q_m values varying with distance. Solid line denotes the fitted regression line. Dashed lines denote upper and lower standard derivation lines of linear regression; (b) Q_m residuals varying with distance. Solid line denotes fitted regression line. Dashed lines denote upper and lower standard derivation lines of linear regression; (c) Q_m values varying with depth; (d) Q_m residuals varying with depth.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Q_m values obtained from ALS station data records at depth range 9~13 km, and its residuals characteristics varying with epicenter distance and depth. (a) Q_m values varying with distance. Solid line denotes the fitted regression line. Dashed lines denote upper and lower standard derivation lines of linear regression; (b) Q_m residuals varying with distance. Solid line denotes the fitted regression line. Dashed lines denote upper and lower standard derivation lines of linear regression; (c) Q_m values varying with depth; (d) Q_m residuals varying with depth.

DownLoad: Full-Size Img PowerPoint

We calculate the averaged Q_m residual by taking the mean value of Q_m residual specified depth range in a one year moving window, which shifts at steps of one month. Figure 9 shows changes of averaged Q_m residual varying with time from SML and ALS stations and specified depth ranges. We can see from the comparison between Figure 9 and Table 1 that significant changes of averaged Q_m residual are highly related to the regional tectonic activity, especially Chi-Chi strong earthquake. Averaged Q_m residual varying with time from SML station located in the source area of the Chi-Chi earthquake for two depth range, 4-9 km and 13-25 km, show obvious tendency rise and lasting higher values fluctuation changes (more than upper standard deviation values) from at least February 1997 until March 1998 and tendency fall fluctuation changes (reach up to lower standard deviation values) from April 1998 until the Chi-Chi earthquake (see Figures 9a and 9b). Averaged Q_m residual varying with time from SML station located in the source area of the Chi-Chi earthquake for depth range more than 25 km show obvious tendency rise and lasting higher values fluctuation changes from January 1998 until the Chi-Chi earthquake (see Figure 9e). Averaged Q_m residual varying with time from ALS station located near the source area of the Chi-Chi earthquake for two depth range, 4-9 km and 13-25 km, show obvious tendency rise and lasting higher values changes from April 1999 until the Chi-Chi earthquake (see Figures 9c and 9d). After the Chi-Chi earthquake, changes of averaged Q_m residual gradually back normal background state. We also see from the comparison between Figure 9 and Table 1 that attenuative dispersion characteristics of direct P-waves in and near the source area of the 17 July 1998 M_W5.6, 17 May 2000 M_W5.6 and 28 July 2000 M_W5.7 (excluded Chi-Chi strong aftershocks in Table 1) have similar patterns to the Chi-Chi earthquake, except that anomalies magnitude have difference.

Figure  9.  Averaged Q_m residuals varying with time from SML and ASL stations and specified depth ranges. The arrows denote time (GMT) of earthquakes occurred more than magnitude M5.0 within epicentral distance of 100 km from SML or ASL station since 1996-2003.

DownLoad: Full-Size Img PowerPoint

Table  1.  M_W ≥ 5.5 earthquakes occurred within a distance of 100 km from SML and ALS observation stations (data from Broadband Array in Taiwan for Seismology)

 | Show Table

DownLoad: CSV

6.   Discussion and Conclusions

We can infer from the results that: ① lasting higher averaged Q_m residual is likely be caused by pre-seismic stress accumulation in and near the source region of Chi-Chi earthquake; ② tendency fall changes of averaged Q_m residual reflects elasticity intensity decrease of source region of Chi-Chi earthquake, which is interpreted as images of a highly fracture and fluid-filled Chi-Chi source region; ③ the difference changes of time period lasting higher averaged Q_m residual from different depth ranges at the same station reflects the difference of shallower and deeper rock properties. The elasticity intensity of deeper medium is higher than that of shallower medium; ④ the difference changes of time period lasting higher averaged Q_m residual from different stations reflects the lateral variation of stress field, the tectonic stress level of source area of Chi-Chi earth- quake is higher than that of farther area from the epicenter of the Chi-Chi earthquake; ⑤The area near observation station with both anomalously increasing and decreasing averaged Q_m residual is likely an unstable environment for future strong earthquake occurrence. In summary, Q_m residual varying with time is a indicator studying component and property changes of the rocks and regional stress field change.

The increase of Q with increasing distance has commonly been observed (e.g., Al-Shukri and Mitchell, 1988; Cong et al, 2000; Liu et al, 2005). Chen et al (2001) reported that P-waves speed was lower in and around the source volume of the Chi-Chi earthquake at a depth range of 5-25 km along the trend of the dipping Chelungpu thrust. Tendency fall changes of averaged Q_m residual from SML station is consistent with lower P-wave speed anomalies in the source area of Chi-Chi earthquake. Li et al (2001) studied the relation between rock porosity permeability property and seismic wave attenuation and propagation velocity and concluded that: ① the velocity was obviously related to porosity, which showed the linear decrease of velocity with increase of porosity, and it was not obviously related to permeability; ② the quality factor was not obviously related to porosity, however it was obviously related to permeability, which showed the exponential decrease of quality factor with increase of permeability. The same studied results about relating P-wave attenuation to permeability was early reported by Akbar et al (1993). Therefore, we consider the fall changes of Q_m is associated with fluid-filled higher density fractures rock volume in the source area of Chi-Chi earthquake, and Q was likely more sensitive than seismic wave speed to the component and property changes of the rocks. Several type of precursors that had been positively identified in relation to the Chi-Chi strong earthquake were found by Liu et al (2004), Chen et al (2004) and Tsai et al (2006).These precursors included the increase of P-wave travel-time residual, the significant ground surface deformation by DInSAR patterns, a significant variation of the b-value and decrease of geomagnetic total field near the Chelungpu fault. These precursors appeared about 2.5 years before the Chi-Chi earthquake. In our study the anomalous changes of averaged Q_m in time and location is quite consistent with the above precursors changes. These different types of precursory phenomena are probably related to accelerating crustal instability just prior to the Chelungpu fault failure, and may be explained by means of the dilatancy and diffusion model for earthquake preparation (Nur, 1972).

We develop a method for measuring intrinsic and attenuative dispersion of direct P-waves. Based on observed group delays, continuous relaxation model and genetic algorithms, we infer Q_m values corresponding to 984 ray paths from SML and ALS stations. The distribution of Q_m are characterized by increase with distance and depth, and is clearly independent of magnitude. Q_m residual from specified depth range which corresponds to one or two layers of the crust block of the Taiwan westward show no obvious correlation with distance, depth and magnitude. The pre-seismic significant changes of averaged Q_m residual about 2.5 years are related to gradual buildup of the east-west tectonic compression stress before the Chi-Chi earthquake in Taiwan. Tendency fall changes of averaged Q_m residual reflects elasticity intensity decrease of regional medium due to the changes of rock porosity and permeability properties. The difference changes of time period lasting higher averaged Q_m residual from different depth ranges at the same station reflects the difference of shallower and deeper rock properties. The elasticity intensity of deeper medium is higher than that of shallower medium. The difference changes of time period lasting higher averaged Q_m residual from different stations reflects the lateral variation of stress field. The tectonic stress level of source area of Chi-Chi earthquake is higher than that of farther area from the epicenter of the Chi-Chi earthquake. This study demonstrate the capability of direct P-waves dispersion for monitoring intrinsic attenuation characteristics and its changes of anelastic medium of the Earth at short distance using waveforms recorded from very small events.