Earthquake Reaearch in China  2019, Vol. 33 Issue (2): 305-319     DOI: 10.19743/j.cnki.0891-4176.201902012
Study of the Change in Wave Velocity Ratio before and after Two Strong Earthquakes Using Airgun Source Data
LIU Zifeng1, YE Beng2, CHEN Jia2, ZHOU Qingyun1, WANG Guangming1, PENG Guanling1, LI Yongli1
1. Yunnan Earthquake Agency, Kunming 650224, China;
2. Office of the Western Yunnan Earthquake Prediction Study Area, China Earthquake Administration, Dali 671000, Yunnan, China
Abstract: Using the signals excited by the large-volume airgun source at the Binchuan transmitting seismic station from January to June, 2016, arrival-time data was acquired at four stations near the epicenter of the Eryuan MS4.5 and MS4.0 earthquakes on February 8, 2016, as well as the epicenter of the Yunlong MS5.0 and Eryuan MS4.6 earthquakes on May 18, 2016 through the waveform cross-correlation technique. The wave velocity ratio of the four stations was calculated using the single-station method. At the same time, the b-value and the focal mechanism consistency parameters of the study area were also calculated. The results show that:(1) the wave velocity ratio of each station experienced a process of decline-recovery-fast rise before the two strong earthquakes, and a significant quasi-synchronous rapid rise occurred within 3-12 days before the earthquake; (2) the timing of the rapid rise of the wave velocity ratio of the four stations before the Yunlong MS5.0 and Eryuan MS4.6 earthquakes were related to the epicentral distance. The station which observed the earliest increase in rapid rise is the farthest one from the epicenter, and the station where the rapid rise appeared in the latest is closest to the epicenter; (3) the form of change of the wave velocity ratio before the earthquake was different between stations located at different directions in the epicenter area; (4) the b-value and the focal mechanism consistency parameter which is commonly used to characterize the stress level both showed a downward trend before the two strong earthquakes, and were consistent with the change in the wave velocity ratio. According to the preliminary analysis, the wave velocity ratio obtained by using airgun source can better reflect the change in the stress state of the underground medium.
Key words: Airgun source     Wave velocity ratio     b-value     Focal mechanism consistency

INTRODUCTION

In the process of evolution of the earth's crust, the physical state of the medium will produce a series of changes, such as micro-fracturing, expansion, plastic hardening and phase transformation. When the seismic wave passes through the crustal medium, the wave velocity will also change accordingly, this is an important basis for studying the physical properties of the media by using wave velocity ratio (Feng Deyi, 1981). Many studies have been conducted in this area by researchers using the Wadati method of multi station and single earthquake, and of single station and multi earthquakes. Feng Deyi et al.(1974, 1978) made a systematic study on the variation characteristics of wave velocity ratio before and after a large number of moderate and strong earthquakes in the 20th century. On the basis of this analysis, many scholars have done a lot of research on the characteristics of wave velocity ratio change before and after significant earthquakes since 2000 (Li Mingxiao et al., 2006; Wang Linying et al., 2008, 2014; Li Yan'e et al., 2014; Weng Zhaoqiang et al., 2015; Liu Zifeng et al., 2018). The results show that before most strong earthquakes or strong aftershock of sequences, the variation process of a normal-decline-recovery-earthquake of the wave velocity ratio was observed, and the wave velocity ratio anomaly generally appeared earlier at the far station and relatively later at the near station. In addition, the similar change pattern was observed before the reservoir-induced earthquakes, that is, the wave velocity ratio also showed a development process of a decline-recovery-earthquake (Zou Zhenxuan et al., 2006; Li Yongli et al., 2012). It can be seen that the existence of wave velocity anomalies before earthquakes has been affirmed by different researchers through studying the characteristics of the wave velocity ratio change before and after strong earthquakes, strong aftershock of sequences and reservoir-induced earthquakes.

The above studies are all aimed at natural earthquakes or reservoir-induced earthquakes. The use of artificial sources to the study of wave velocity change is mainly to indirectly characterize the change of subsurface wave velocity through the relative change of travel time about specific seismic phases. With the development of travel time measurement technology, some progress has been made in the indirect and accurate measurement of wave velocity changes (Silver P. G. et al., 2007; Wang Baoshan et al., 2008; Wang Weitao et al., 2009; Xu Hui et al., 2015; Wei Yunyun et al., 2016; Chen Jia et al., 2017; Zhang Yuansheng et al., 2017; Zhou Qingyun et al., 2018; Liu Xuzhou, 2018). The results show that there are indeed changes in the travel time of specific phases before many strong earthquakes, and it is inferred that the change in travel time is caused by the change of wave velocity, and the wave velocity variation is mainly caused by the change in the physical state of the subsurface medium with the stress state.

The airgun source has rich phases, and clear P-wave and S-wave phases can be obtained by stacking the records of multiple shots (Chen Yong et al., 2017; Zhang Yuansheng et al., 2017). By utilizing the arrival time data of P-wave and S-wave phases of the airgun source, the wave velocity ratio can then be directly calculated using the single station method. Therefore, based on the airgun signatures excited by the large-volume airgun source of the Binchuan transmitting seismic station, this paper calculates the wave velocity ratio of the stations near the epicenters before and after the two strong earthquakes in the Yunlong-Eryuan area in 2016. In order to further analyze whether the change of the wave velocity ratio calculated using the airgun source signatures can truly reflect the micro-dynamic change of physical properties of the underground medium. The paper also calculates the b-value of the seismological parameter, which is commonly used to characterize the stress level (Wyss M. et al., 2000; Yi Guixi et al., 2008, 2010, 2013), and the focal mechanism consistency (Chen Yong, 1978; Michael A. J. et al., 1990; Lu Zhong et al., 1997) to corroborate each other.

1 DATA SELECTION AND PRE-PROCESSING 1.1 The Active Seismic Source Data of Binchuan Transmitting Seismic Station

The large-volume airgun signals transmitting seismic station in Binchuan, Yunnan was built in 2011 and officially put into operation in 2012. The signal is continuously excited within the constraints of water level and airgun status. This paper selects the airgun excitation signals of the transmitting seismic station from January to June 2016. Fig. 1 shows the daily frequency of the airgun excitation during the study period. It can be seen that the excitation can be guaranteed about 13 times a day and the data is basically continuous in this period. This provides relatively continuous stable data for analysis of the changes in wave velocity ratio before the two strong earthquakes in January and June, 2016. In order to study the characteristics of the wave velocity ratio change before the strong earthquakes, the airgun signals recorded by the nearest four stations to the earthquakes are selected for calculation and analysis. The spatial distribution of the selected stations, epicenters and airgun are shown in Fig. 2. The four stations well encircle the epicenters and locate respectively in the four quadrants of the epicenter (Fig. 2), which provides a convenient condition for a comparative analysis of wave velocity ratio variation in different azimuths in the epicenter area.

 Fig. 1 Daily frequency of airgun excitation at Binchuan transmitting seismic station January 1-June 30, 2016

 Fig. 2 Distribution of epicenters of MS≥4.0 earthquakes and the seismic stations selected in this study

In this paper, the airgun source signals received by the four stations are used to pick up the arrival time of P and S waves of the corresponding station. The airgun signal waveform received by each station carries the interference information such as zero drift of instrument and unstableness of the base. Therefore, it is necessary to perform de-meaning, de-trending and anti-tilting to the airgun signals before picking up the phase data. Considering that the dominant frequency of airgun is concentrated at 2-5Hz (Chen Jia et al., 2016), the waveform is pre-processed by a 2-5Hz bandpass filter, and finally, the waveform data with significant interference is removed manually.

In this paper, the airgun excitation moment is the origin time of the earthquake, but the accuracy of the artificially recorded excitation moment is only 1s. The low accuracy will affect the calculation result. Therefore, the P-wave arrival time of ckt is used as the excitation moment (Fig. 3(a)), for the distance of ckt is only 50m to the airgun source which can be approximated as source information (Chen Jia et al., 2017). The P-wave and S-wave arrival-times are intercepted respectively in the vertical and horizontal direction with a certain window length as a template, and the excitation records are superimposed and then cross-correlated with the waveform each day, and finally, according to the correlation coefficient, the arrival time of the corresponding phase is determined. This study obtains the arrival time data of the 4 stations with 0.8 as the lower limit of the correlation coefficient. The number of records in which both P-wave and S-wave arrival time can be picked up is 122 in station No.52536, 115 in station YUL, 98 in station EYA, and 76 in station EY16.

 Fig. 3 Sketch diagram of the excitation moment(a), and P-wave and S-wave arrival times(b)
1.2 Earthquake Catalogue

This paper uses the earthquake catalogue from January to June 2016 provided by the Yunnan Digital Seismic Network to calculate the b-value of Yunlong-Eryuan and its surrounding areas. In order to avoid the influence of data selection on the statistical results, the removal of the aftershock sequence, de-clustering of earthquake swarms and the analysis of completeness of the catalogue are crucial. This paper manually performed a second check based on the preliminary de-clustering with the kk method, and finally obtained 677 earthquakes of magnitude 0 or higher. For the above catalogue, the EMR method provided by Zmap software (Wiemer S., 2001) is used to calculate the minimum magnitude of completeness, MC, of the study area during the study period, which is MC=1.0. The epicenter distribution of the earthquake with ML≥1.0 in the study area from January to June 2016 is shown in Fig. 4.

 Fig. 4 Epicenter distribution(ML≥1.0)
1.3 Focal Mechanism Solutions

The focal mechanism solutions used to calculate the consistency parameter are integrated by the relevant results obtained by the CAP method and the P-wave initial motion method. In order to eliminate the influence of sequence aftershocks on the consistency parameters of the focal mechanism, this paper removes the sequence aftershocks and finally obtains the focal mechanism solutions for 57 ML≥2.5 earthquakes in Yunlong-Eryuan and its adjacent areas in 2015 and 2016 (Fig. 5), including 32 events with magnitude 2.5 to 2.9, 17 with magnitude 3.0 to 3.9, 7 with magnitude 4. 0 to 4.9, and one with magnitude 5.0 to 5.9.

 Fig. 5 Focal mechanism distribution (ML≥2.5)
2 PRINCIPLE OF THE CALCULATION METHOD 2.1 Wave Velocity Ratio Calculation Method

In this study, the simplified method proposed by Li Yongli et al. (2012) is used to calculate the wave velocity ratio. This method assumes that the distance from the seismic source to station is Δ (in this paper, Δ is the distance between the airgun and station), and the P-wave velocity is vP, the travel time is tP, S-wave velocity is vS, the travel time is tS. Then the wave velocity ratio is:

 $\frac{{{v_{\rm{P}}}}}{{{v_{\rm{S}}}}} = \frac{{\Delta /{t_{\rm{P}}}}}{{\Delta /{t_{\rm{S}}}}} = \frac{{{t_{\rm{S}}}}}{{{t_{\rm{P}}}}}$ (1)

where, tS=TS-T0, tp=Tp-T0; Tp, TS are respectively the P-wave and S-wave arrival time, T0 is the origin time of earthquake(in this paper, T0 is the excitation moment of airgun), Substituting them into equation (1) and by simplification, we get:

 $\frac{{{v_{\rm{P}}}}}{{{v_{\rm{S}}}}} = 1 + \frac{{{T_{\rm{S}}} - {T_{\rm{P}}}}}{{{T_{\rm{P}}} - {T_0}}}$ (2)

Using equation(2), the wave velocity ratio, vP/vS, can be calculated.

The paper takes one excitation recorded at station EYA as an example to illustrate the calculation process. According to ckt records, an excitation was conducted at 21:59:59.780 on January 7, 2016 (namely T0), as shown in the records of station EYA, the P-wave arrival time of this excitation is at 22:00:11. 983(namely Tp), the S-wave arrival time is at 22:00:20.401(namely TS), substituting T0, TP and TS into equation(2) yields a wave velocity ratio of 1.6898. According to the above method, we can obtain the wave velocity ratio of each excitation for each station, and on this basis, analyze its time-series variation characteristics. In order to reduce the effect of random variations on the calculation results, the Savitzky-Golay smoother integrated in Origin Software was used to smooth the calculation results of wave velocity ratio of each station in the time-series variation analysis.

The Savitzky-Golay smoother was originally proposed by Savitzky A. and Golay M. in 1964 (Cai Tianjing et al., 2011) and is widely used for data stream smoothing and denoising. It is based on a polynomial in the time domain and is solved by moving the window using least squares method. It is the best fitting method, which can well preserve the distribution characteristics of the relative maxima, minima and width of the data stream, and can ensure that the original data is not distorted to the utmost. The specific process of the Savitzky-Golay smoothing method is as follows:

Let a set of data in x(n) be x(i), i=-M, …., 0…., M, by constructing a p-order polynomial fi to fit x(i), where,

 ${f_i} = {a_0} + {a_1}i + {a_2}{i^2} + \cdots + {a_p}{i^p} = \sum\limits_{k = 0}^p {{a_k}} {i^k}p \le 2M$ (3)

and its total squared error and E are:

 $E = \sum\limits_{i = - M}^M {{{\left[ {{f_i} - x(i)} \right]}^2}} = \sum\limits_{i = - M}^M {{{\left[ {\sum\limits_k^p {{a_k}} {i^k} - x(i)} \right]}^2}}$

In order to make the curve have the highest degree of fit and the minimum fitting error, let the derivative of E to each coefficient be 0, so as to achieve the total squared error and minimum E, namely,

 $\frac{{\partial E}}{{\partial {a_r}}} = 0, r = 0, 1, 2, \cdots, p$

Hence, we get:

 $\frac{{\partial E}}{{\partial {a_r}}} = \frac{\partial }{{\partial {a_r}}}\left\{ {\sum\limits_{i = - M}^M {{{\left[ {\sum\limits_{k = 0}^p {{a_k}} {i^k} - x(i)} \right]}^2}} } \right\} = 2\sum\limits_{i = - M}^M {\left[ {\sum\limits_{k = 0}^p {{a_k}} {i^k} - x(i)} \right]} {i^r} = 0$ (4)
 $\sum\limits_{k = 0}^p {{a_k}} \sum\limits_{i = - M}^M {{i^{k + r}}} = \sum\limits_{i = - M}^M x (i){i^r}$ (5)

let ${F_r} = \sum\limits_{i = - M}^M x (i){i^r}, {S_{k + r}} = \sum\limits_{i = - M}^M {{i^{k + r}}}$, then equation (5) can be written as:

 ${F_r} = \sum\limits_{k = 0}^p {{a_k}} {S_{k + r}}$ (6)

Therefore, by selecting M (the number of single-sided points to be fitted) and p (polynomial order), Fr and Sk+r can be obtained, and substituting into equation (6), the coefficients a0, a1, …. ap can be obtained. Thus, the polynomial fi is determined, and the smoothed data value is obtained.

In this paper, the Savitzky-Golay smoother is used to smooth the wave velocity ratio of each station obtained by the single-station method. Fig. 6 shows the fitting smooth curve of the station EYA with the single side points M=10 and the polynomial order p=2. It can be seen from Fig. 6 that the smoothing can effectively remove the random step change, which is convenient to analyze the overall change trend of the wave velocity ratio.

 Fig. 6 Comparison of the wave velocity ratio of station EYA before and after smoothing
2.2 A Brief Introduction of the Methods for Calculating b-value

Gutenberg-Ricker's law proposes that the magnitude is linear with the logarithm of the corresponding frequency, i.e. lgN=a+bM, where N represents the frequency of earthquakes above magnitude M and the a-value represents the level of seismic activity in the statistical region, and the b-value reflects both the proportional relationship between small earthquakes and large earthquakes, and the stress state of the medium, explaining why it is widely used in seismic risk analysis. Generally, the decrease of the b-value is considered to be the expression of the enhancement of regional stress and seismic risk (Li Yongli et al., 2002; Yi Guixi et al., 2008, 2010, 2013; Wang Hui et al., 2012; Mao Yan et al., 2016). The most commonly used methods for calculating the b-value are the least squares method and the maximum likelihood method. This study uses the maximum likelihood method, which is the default method provided by Zmap package (Wiemer S., 2001).

The maximum likelihood method is a commonly used point estimation method. Its essence is to take the parameter value that maximizes the probability of occurrence of the observed earthquake sample as an estimate of the unknown parameter (Chen Yang et al., 2013). The calculation formula (Aki K., 1965) is expressed as:

 $\hat b = \frac{{{\rm{ }}N{\rm{lge }}}}{{\sum\limits_{i = 1}^N {\left({{M_i} - {M_0}} \right)} }}$ (7)

Where, $\hat{b}$ represents the estimate of b-value, N the total numbers of earthquake, Mi the magnitude of the ith earthquake, and M0 the threshold magnitude.

2.3 A Brief Introduction of the Methods for Calculating Focal Mechanism Consistency Parameters

The use of a large number of focal mechanism solution data can infer the characteristics of the regional stress field. Seismologists have proposed many methods in this regard (Michael A. J., 1984, 1987; Gephart J. W. et al., 1984; Xu Zhonghuai et al., 1987; Zhong Jimao et al., 2006). Although the occurrence of small and moderate earthquakes is highly random and it is inconvenient to conduct specific analysis individually, Hardebeck et al. (2006) pointed out that the direction of the stress tensor can be accurately constrained according to a large number of discretely distributed data. To solve the characteristics of such spatially non-uniform stress fields in the focal mechanism solutions, Michael proposed a superposition stress field inversion method (SSI), which simulates the distribution of non-uniform stress fields by superimposing disturbances on a uniform stress field to obtain the spatial distribution and relative size φ=(S2-S3)/(S1-S3) as well as the inversion variance of the three principal stress axes S1, S2 and S3 (i.e. the principal compressional stress axis, the intermediate principal stress axis and the extensional principal stress axis, respectively). The inversion variance is defined as the mean of the misfit angle (the angle between the sliding vector of a single earthquake and the theoretical sliding vector produced under the assumed stress tensor) and the sum of the squares of the difference between them and the mean. It is a quantitative indicator to measure the degree of consistency between the stress field released from the seismic source and the regional tectonic stress field (Michael A. J., 1987; Michael A. J. et al., 1990). When the variance is < 0.1, it means that a unified stress tensor can be used to explain the observed focal mechanism solution. It can also be understood that the stress field in the region is uniform and the focal mechanism is consistent. When the variance is >0.2, it indicates that the stress field in this region is heterogeneous in time and space, or the focal mechanism of the region is relatively disordered (Lu Zhong et al., 1997). Experiments have shown that the regional stress tensor obtained by the above method is in line with the observation facts. Many scholars at home and abroad have carried out a lot of research based on this method (Michael A. J., 1991; Wiemer S. et al., 2002; Zhang Zhiwei et al., 2015; Li Jin et al., 2015). Therefore, this paper uses the Michel method integrated by the Zmap software package (Weimer S., 2001) to calculate the focal mechanism consistency parameters of the study area.

3 ANALYSIS OF THE CALCULATION RESULTS 3.1 Analysis of the Variation Characteristics of Wave Velocity Ratio before and after Moderate-strong Earthquakes

The wave velocity ratio before the MS4.5 and MS4.0 earthquakes on February 8, 2016 is calculated at three stations, namely station No.53036, EYA and YUL. Station EY16 was not yet set up then. Fig. 7 shows that before the earthquakes, the wave velocity ratio of the 3 stations all showed a process of decline-recovery, but the range and timing of the change were different. At station 53036, the wave velocity ratio began to recover rapidly on February 5 from a continuous declining trend; at station EYA, the wave velocity ratio was in a fluctuant downward trend during the period from January 4 to February 1, and the rate of decline began to increase on January 18, and the turning point appeared on January 22, with the recovery rate speeding up on February 5. At station YUL, the wave velocity ratio began to turn upward after a fluctuant decline from January 4 to January 28. It is not difficult to find that the decline of the wave velocity ratio of stations EYA and YUL before the Eryuan earthquake on February 8 was significantly higher than that of station No.53036. The increased range of station YUL immediately before the earthquake was more significant than that of station 53036 and EYA. During the change of the wave velocity ratio, the rapid decline appeared synchronously at stations YUL and EYA, while the increased recovery appeared almost synchronously at stations 53036 and EYA before the earthquake, and the change patterns were very similar. After the MS4.5 and MS4.0 earthquakes in Eryuan on February 8, 2016, the wave velocity ratio of stations 53036, EYA and YUL began to decrease quasi-synchronously.

 Fig. 7 Time series curve of wave velocity ratio variation of each station

Due to equipment failure, there was a lack of readings from February 23 to March 17, so, the analysis was focused mainly on the variation characteristics of wave velocity ratio around the Yunlong MS5.0 and Eryuan MS4.6 earthquakes on May 18, 2016. Since March 18, the wave velocity ratio of station 53036 experienced a change process of small-scale slow decline - rebound - accelerated recovery - rapid decline - earthquake; the change of the ratio in station YUL shows the characteristics of trending decline - slow recovery - turning and slow decline - turning and rapid rise - turning and rapid decline-earthquake; while the change in the wave velocity ratio of station EY16 experienced a process of slow decline - turning and rapid rise - earthquake; among the 4 stations, the wave velocity ratio change in station EYA is the more intensive. After two rapid rise and decline processes, the ratio turned to rise on April 20, and on May 6, the rising rate began to speed up. It is clear from Fig. 7 that the wave velocity ratio evolution pattern of decline-turning to rise the rising rate increasing immediately before earthquake is the common variation of the four stations before the Yunlong earthquake, which is consistent with the trend of the pre-earthquake wave velocity ratio obtained by using natural earthquakes (Li Mingxiao et al., 2006; Wang Linying et al., 2008, 2014; Li Yan'e et al., 2014; Weng Zhaoqiang et al., 2015; Liu Zifeng et al., 2018). Fig. 7 shows that the rise of the wave velocity ratio appeared earliest in station 53036, followed by station EYA. The anomaly observed by station YUL appeared slightly later than station EYA, and the station EY16 is the latest to observe the rise. It can be seen that the rapid rise of the wave velocity ratios of the four stations has a certain relation with the epicentral distance. Station 53036 where the increase of rising rate appeared earliest is the farthest station to the epicenter (40km), while station EY16 which is the nearest to the epicenter (12km) is the latest to observe the rapid rise, which is consistent with the conclusions obtained by Weng Zhaoqiang et al. (2015) and Liu Zifeng et al. (2018) using natural earthquakes. After the earthquakes on Yunlong MS5.0 and Eryuan MS4.6 on May 18, 2016, the wave velocity ratio of the four stations fluctuated greatly, and the variation trend was different, which may be related to the post-earthquake stress adjustment.

Comparing the characteristics of the wave velocity ratios of the stations in different azimuths in the epicenter area before the MS4.5 and MS4.0 earthquakes on February 8, 2016, it is found that in the rapid decline process, the change pattern of wave velocity ratio is basically the same in the stations EYA and YUL located respectively at the first and third quadrant, while in the rapid rise process, the change pattern of the stations EYA and 53036 located respectively on the 2nd and 4th quadrant are very much the same. The change pattern of wave velocity ratios is relatively consistent in stations EY16 and 53036 located respectively in the 2nd and 4th quadrants. The change in station EY16 that passes through the epicenter is more notable than that of station 53036 without passing the epicenter. The changes in stations EYA and YUL in the 1st and 3rd quadrant are basically the same, and the change in station EYA is larger than that of the station YUL.

The wave velocity ratios before the two strong earthquakes experienced a process of decline-recovery-recovery speed-up, but there were some differences in the time and range of the change. The DD mode mechanism of "expansion-water inflowing" can reasonably explain the development process of the decline-rise-earthquake of the pre-earthquake wave velocity, that is, the wave velocity decreases with the increase of fracture density in the dry fracture environment, and increases in a water-saturated fracture setting (Feng Deyi, 1981). According to this analysis, the decrease of the wave velocity ratio before the two strong earthquakes may be caused by the increase of fracture density caused by stress enhancement during the earthquake preparation. The rebound process before the earthquake may result from the fluid infiltration in the late stage of the earthquake preparation, and the participation of the fluid leads to the cracks reaching a water saturation state, which increases the wave velocity and the wave velocity ratio, vP/vS, as well.

3.2 Analysis of the Characteristics of b-value as a Function of Time

Using the maximum likelihood method and the official catalogue provided by the Yunnan Digital Seismic Network with a threshold magnitude of ML1.2 (>MC), the daily sliding b-value of accumulatively 15 days from January to June 2016 of the Yunlong-Eryuan and adjacent areas is calculated under the condition that a sample size of 15 is met. The calculation results are shown in Fig. 8. The b-value before the two strong earthquakes dropped by a big margin, and after the earthquakes, it started to rise after a short period with a low b-value. According to the theory that the stress is inversely proportional to the b-value (Scholz C. H., 1968; Wyss M. et al., 2000), it can be inferred that the decrease in the b-value before the moderate-strong earthquake indicates the increase of stress, and the stress level decreases with the release of energy after an earthquake, so the b-value began to rise after the earthquakes.

 Fig. 8 Time series curve of b-value of Yunlong-Eryuan and its adjacent areas
3.3 Analysis of Variation Characteristics of Focal Mechanism Consistency Parameter with Time

For the focal mechanism solutions of 57 ML≥2.5 earthquakes in Yunlong-Eryuan and its adjacent areas in 2015 and 2016, we use the Michael method and take 8 earthquakes as the window length and 5 earthquakes as the step length to calculate the variations of stress tensor variance over time. Fig. 9 shows the variation of the variance over time from January to June 2016. It can be seen that the variance value before the two strong earthquakes showed a downward trend. Michael A. J. et al., (1990) suggested that the size of the misfit angle or stress tensor variance can characterize the consistency degree of the focal mechanism. The low value reflects a good consistency of the focal mechanism, and a uniform stress field means high stress levels; while the high value reflects a poor consistency of the focal mechanism and the stress field is non-uniform. According to the analysis, the value of the Variance of the focal mechanism consistency parameter before the two strong earthquakes showed a downward trend, which is a manifestation of stress enhancement in the source region.

 Fig. 9 Time series curve of the value of Variance of the focal mechanism consistency parameter in Yunlong-Eryuan and its adjacent areas
4 ANALYSIS ON THE RELIABILITY AND STABILITY OF THE CALCULATION RESULTS

In calculating the wave-velocity ratio using the single-station method, the main factors affecting the calculation results are the pick-up accuracy of arrival time of P-wave and S-wave, and the accuracy of the earthquake location. In this study, the airgun source is used to calculate the wave velocity ratio, and the source location is fixed. Therefore, the influence of earthquake location accuracy on the calculation results can be ruled out. The main influencing factor is the pick-up accuracy of the arrival time of direct P-wave and S-wave. In this paper, the waveform data is pre-processed before the arrival time is picked up by the waveform cross-correlation method. Finally, with high-quality waveform data, the arrival times are picked up taking 0.8 as the lower limit of the correlation coefficient. Accuracy can reach 10-2-10- 3s, the absolute arrival time difference error obtained is small, and the average standard deviation of the absolute arrival-time difference of P-wave and S-wave is 0.14 and 0.19, respectively. In addition, in analyzing the temporal variation characteristics of vP/vS, in order to reduce the influence of random variation on the calculation results, this paper adopts the wave velocity ratios excited by 10 excitations for the moving average processing, which can effectively improve the stability and reliability of vP/vS over time. Finally, the standard deviation of the wave velocity ratio is 0.0076, 0.0072, 0.0173 and 0.0106 for stations EY16, YUL, EYA and 53036, respectively, showing a good stability of the calculated wave velocity ratios.

5 CONCLUSION AND DISCUSSION

Based on the excitation signals of the large-capacity airgun source in the Binchuan transmitting seismic station from January to June 2016, this paper acquired the absolute arrival time data of P-wave and S-wave of the four stations near the epicenters of the Eryuan MS4.5 and MS4.0 earthquakes on February 8, 2016 and the Yunlong MS5.0 and Eryuan MS4.6 earthquakes on May 18, 2016, using the waveform cross-correlation method. The wave-velocity ratios of the four stations are calculated by the single-station method, and at the same time, the b-value, the focal mechanism consistency parameter, and variance of Yunlong-Eryuan and its adjacent areas in the same periods are obtained. The main conclusions are as follows:

(1) The wave velocity ratios of the four stations before both strong earthquake groups experienced a process of decline-recovery-fast rise, and the time and range of the change were different. The quasi-synchronization of rapid rise occurring within 3-12 days before the earthquakes was prominent.

(2) The time of the rapid rise of the wave velocity ratio appearing in the four stations before the Yunlong earthquake is related to the epicentral distance to a certain extent. The rising rate increased earliest in station 53036 which is the farthest from the epicenter, and station EY16, which is the nearest the epicenter, is the latest to observe the rapid rise, and similar phenomena has been found in related studies conducted using natural earthquakes.

(3) Comparing the characteristics of wave velocity changes of the four stations located in different azimuths in the epicenter area before the Yunlong earthquake on May 18, 2016, it is found that the change pattern in station EY16 and 53036 located in the 2nd and 4th quadrants are consistent, and the change in station EY16 which runs through the epicenter was more striking than that of station 53036 which does not run through the epicenter. The change patterns of station EYA and YUL located in the 1st and 3rd quadrants are basically the same, and the change amplitude in station EYA was larger than that of station YUL. However, before the MS4.5 and MS4.0 earthquakes in Eryuan on February 8, 2016, the change pattern of the wave velocity ratio in a declining process was basically the same in the stations EYA and YUL which are located in the 1st and 3rd quadrants, and in the rapid rise process, the change pattern of stations EYA and 53036 located at the 1st and 4th quadrant was very much similar. It can therefore be seen that there are differences in the wave velocity ratio changes of stations located in different azimuths, but the changes did not show obvious regularity.

(4) The b-value and the focal mechanism consistency parameter, variance, commonly used to characterize the stress level were calculated. The b-value and variance value of the study area before strong earthquakes were significantly decreased, which is a manifestation of significant increase in stress levels near the epicenter.

The comprehensive analysis shows that the variation characteristics of the wave velocity ratio before strong earthquakes obtained using the airgun source are consistent with the relevant research results using natural earthquakes, and the trends and patterns of the wave velocity ratio, b-value and Variance value around the two strong earthquakes are consistent. This result indicates that the change in the wave velocity ratio calculated using the airgun source can be used to characterize the stress state and provide some criteria for earthquake situation tracking.

ACKNOWLEDGEMENT

The author would like to express their heartfelt thanks to all the institutions and individuals who have contributed to the large-volume airgun experiment in Binchuan.

REFERENCES
 Aki K. Maximum likelihood estimate of b in the formula logN=a-bM and its confidence limits[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 1965, 43(2): 237-239. Cai Tianjing, Tang Han. Review on the least-squares-fit principle of Savitzky-Golay smoothing filter[J]. Digital Communication, 2011, 38(1): 63-68, 82 (in Chinese with English abstract). Chen Jia, Li Xiaobin, Yang Jun, Ye Beng. Large volume air-gun source spectrum characteristics of Binchuan, Yunnan[J]. Earthquake Research in China, 2016, 32(2): 216-221 (in Chinese with English abstract). Chen Jia, Ye Beng, Gao Qiong, Wang Jun, Li Xiaobin, Yang Jun. Research on variation characteristics of wave velocity ratio before and after the 2016 Yunlong MS5.0 earthquake using airgun signals[J]. Journal of Seismological Research, 2017, 40(4): 550-556, 677 (in Chinese with English abstract). Chen Yang, Lv Yuejun, Xie Zhuojuan, Pan Long. Review of the study of seismicity parameter b-value[J]. Bulletin of the Institute of Crustal Dynamics, 2013, 38-47 (in Chinese with English abstract). Chen Yong. Consistency of focal mechanism as a new parameter in describing seismic activity[J]. Chinese Journal of Geophysics, 1978, 21(2): 142-159 (in Chinese with English abstract). Chen Yong, Wang Baoshan, Yao Huajian. Seismic airgun exploration of continental crust structures[J]. Science China Earth Science, 2017, 60(10): 1739-1751. DOI:10.1007/s11430-016-9096-6 Feng Deyi, Tan Aina, Wang Kefen. Velocity anomalies of seismic waves from near earthquakes and earthquake prediction[J]. Acta Geophysica Sinica, 1974, 17(2): 84-98 (in Chinese with English abstract). Feng Deyi, Gu Jinping, Wang Zhouyuan, Sheng Guoying, Luo Ruiming, Li Kelong. Preliminary study of the velocity anomalies of seismic waves before and after some strong and moderate earthquakes in western China (Ⅲ):Variations of curvature of the Wadati diagrams[J]. Acta Geophysica Sinica, 1978, 21(4): 292-309 (in Chinese with English abstract). Feng Deyi. Velocity Anomaly of Seismic Waves[M]. Beijing: Seismological Press, 1981. 2-15, 140-149, 236-246 (in Chinese). Gephart J.W., Forsyth D.W. An improved method for determining the regional stress tensor using earthquake focal mechanism data:Application to the San Fernando earthquake sequence[J]. Journal of Geophysical Research:Solid Earth, 1984, 89(B11): 9305-9320. DOI:10.1029/JB089iB11p09305 Hardebeck J.L., Michael A.J. Damped regional-scale stress inversions:Methodology and examples for southern California and the Coalinga aftershock sequence[J]. Journal of Geophysical Research:Solid Earth, 2006, 111(B11): B11310. Li Jin, Zhou Longquan, Long Haiying, Nie Xiaohong, Guo Yin. Spatial-temporal characteristics of the focal mechanism consistency parameter in Tianshan (within Chinese territory) seismic zone[J]. Seismology and Geology, 2015, 37(3): 792-803 (in Chinese with English abstract). Li Mingxiao, Liu Jie. Study on velocity ratio (vP/vS) anomaly of earthquake sequences in Yunnan region[J]. Earthquake, 2006, 26(1): 26-34 (in Chinese with English abstract). Li Yan'e, Wang Linying, Zheng Xuyao. Restudy of the variation of vP/vS before and after the Wenchuan earthquake[J]. Acta Seismologica Sinica, 2014, 36(3): 425-432 (in Chinese with English abstract). Li Yongli, Cai Jingguan, Cao Ke. Modulating ratio and b value in strong seismicity course in Yunnan region[J]. Journal of Seismological Research, 2002, 25(1): 25-30 (in Chinese with English abstract). Li Yongli, Zhao Xiaoyan, Fu Hong. Research on variation of wave velocity ratio before and after Water Storage for Xiaowan Reservoir[J]. Journal of Seismological Research, 2012, 35(4): 464-470 (in Chinese with English abstract). Liu Xuzhou. Study on Seismic Travel-Time Variation of Qilianshan Region Using Active Source Data[D]. Lanzhou Institute of Seismology, China Earthquake Administration, 2018 (in Chinese with English abstract). Liu Zifeng, Zhang Tianji, Fu Hong, Li Zhirong, Wang Guangming, Zhao Xiaoyan, Yu Dayuan. Characteristics of wave velocity ratio change before Yunnan Tonghai MS5.0 earthquake in 2018[J]. Journal of Seismological Research, 2018, 41(4): 494-502 (in Chinese with English abstract). Lu Zhong, Wyss M., Pulpan H. Details of stress directions in the Alaska subduction zone from fault plane solutions[J]. Journal of Geophysical Research:Solid Earth, 1997, 12(B3): 5383-5402. Mao Yan, Liu Zifeng, Ye Jianqing, Li Zhonghua. Analysis on strong earthquake risk of Xiaojiang fault zone[J]. Journal of Seismological Research, 2016, 39(2): 213-217 (in Chinese with English abstract). Michael A.J. Determination of stress from slip data:Faults and folds[J]. Journal of Geophysical Research:Solid Earth, 1984, 89(B13): 11517-11526. DOI:10.1029/JB089iB13p11517 Michael A.J. Use of focal mechanisms to determine stress:A control study[J]. Journal of Geophysical Research:Solid Earth, 1987, 92(B1): 357-368. DOI:10.1029/JB092iB01p00357 Michael A.J., Ellsworth W.L., Oppenheimer D.H. Coseismic stress changes induced by the 1989 Loma Prieta, California earthquake[J]. Geophysical Research Letters, 1990, 17(9): 1441-1444. DOI:10.1029/GL017i009p01441 Michael A.J. Spatial variations in stress within the 1987 Whittier Narrows, California, aftershock sequence:New techniques and results[J]. Journal of Geophysical Research:Solid Earth, 1991, 96(B4): 6303-6319. DOI:10.1029/91JB00195 Scholz C.H. The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[J]. Bulletin of the Seismological Society of America, 1968, 58(1): 399-415. Silver P.G., Daley T.M., Niu F.L., Majer E. L. Active source monitoring of cross-well seismic travel time for stress-induced changes[J]. Bulletin of the Seismological Society of America, 2007, 97(1B): 281-293. DOI:10.1785/0120060120 Wang Baoshan, Zhu Ping, Chen Yong, N iu, Fenglin, Wang Bin. Continuous Subsurface velocity measurement with coda wave interferometry[J]. Journal of Geophysical Research:Solid Earth, 2008, 113(B12): B12313. DOI:10.1029/2007JB005023 Wang Hui, Cao Jianling, Jing Yan, Li Zhen. Spatiotemporal pattern of b-value before major earthquakes in the Sichuan-Yunnan region[J]. Seismology and Geology, 2012, 34(3): 531-543 (in Chinese with English abstract). Wang Linying, Guo Yongxia, Liu Fang, Jiang Changsheng. Temporal vP/vS variation characteristics in different zones of China's Capital Circle area before and after Wen'an earthquake[J]. Acta Seismologica Sinica, 2008, 21(3): 243-257 (in Chinese with English abstract). DOI:10.1007/s11589-008-0243-z Wang Linying, Li Yan'e, Zheng Xuyao, Wang Shengwen. Temporal variation of vP/vS at single seismic station before the 2013 Lushan MS7.0 earthquake[J]. Acta Seismologica Sinica, 2014, 36(1): 42-58 (in Chinese with English abstract). Wang Weitao, Wang Baoshan, Ge Hongkui, Chen Yong, Yuan Songyong, Yang Wei, Li Yijin. Using active source to monitor velocity variation in shallow sediment caused by the Wenchuan earthquake[J]. Earthquake Research in China, 2009, 25(3): 223-233 (in Chinese with English abstract). Wei Yunyun, Wang Haitao, Su Jinbo, Chen Xiangjun, Wang Qiong. The preliminary study on travel time abnormal variation of reflection wave phase of air-gun in Xinjiang before two earthquakes with MS5.0[J]. Earthquake Research in China, 2016, 32(2): 270-281 (in Chinese with English abstract). Weng Zhaoqiang, Liang Xiangjun, Wu Haoyu, Liu Linfei, Li Li. The variation of vP/vS at single seismic station of Shanxi region studied by Wadati method of single station and multi-earthquake[J]. Earthquake Research in China, 2015, 31(1): 89-100 (in Chinese with English abstract). Wiemer S. A software package to analyze seismicity:ZMAP[J]. Seismological Research Letters, 2001, 72(2): 373-382. Wiemer S., Gerstenberger M., Hauksson E. Properties of the aftershock sequence of the 1999 MW7.1 Hector Mine earthquake:Implications for aftershock hazard[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 1227-1240. DOI:10.1785/0120000914 Wyss M., Schorlemmer D., Wiemer S. Mapping asperities by minima of local recurrence time:San Jacinto-Elsinore fault zones[J]. Journal of Geophysical Research:Solid Earth, 2000, 105(B4): 7829-7844. DOI:10.1029/1999JB900347 Xu Hui, Liu Xuejun, Wang Bin, Wang Baoshan. Monitoring temporal variation of subsurface seismic velocity in Xiaojiang fault zone using direct wave cross correlation time delay technique with active source[J]. Journal of Seismological Research, 2015, 38(1): 7-15 (in Chinese with English abstract). Xu Zhonghuai, Wang Suyun, Huang Yurui, Gao Ajia, Jin Xiaofeng, Chang Xiangdong. Directions of mean stress axes in southwestern China deduced from microearthquake data[J]. Acta Geophysica Sinica, 1987, 30(5): 476-486 (in Chinese with English abstract). Yi Guixi, Wen Xueze, Su Youjin. Study on the potential strong-earthquake risk for the eastern boundary of the Sichuan-Yunnan active faulted-block, China[J]. Chinese Journal of Geophysics, 2008, 51(6): 1719-1725 (in Chinese with English abstract). Yi Guixi, Wen Xueze, Zhang Zhiwei, Long Feng, Ruan Xiang, Du Fang. Study on potential strong earthquake risk in Mabian area, southern Sichuan[J]. Seismology and Geology, 2010, 32(2): 282-293 (in Chinese with English abstract). Yi Guixi, Wen Xueze, Xin Hua, Qiao Huizhen, Wang Siwei, Gong Yue. Stress state and major-earthquake risk on the southern segment of the Longmen Shan fault zone[J]. Chinese Journal of Geophysics, 2013, 56(4): 1112-1120 (in Chinese with English abstract). Zhang Yuansheng, Wang Baoshan, Chen Yong, Wang Lanmin, Zou Rui, Qin Manzhong, Guo Xiao, Shen Xuzhang, Wei Congxin, Liu Xuzhou, Wang Yahong, Sun Dianfeng, Guo Yingxia, Yin Liang. Travel-time variations before and after two major earthquakes derived from active-source seismic data[J]. Chinese Journal of Geophysics, 2017, 60(10): 3815-3822 (in Chinese with English abstract). Zhang Zhiwei, Zhou Longquan, Long Feng, Ruan Xiang. Spatial and temporal characteristic of stress field for Wenchuan MS8.0 and Lushan MS7.0 earthquake sequence[J]. Seismology and Geology, 2015, 37(3): 804-817 (in Chinese with English abstract). Zhong Jimao, Cheng Wanzheng. Determination of directions of the mean stress field in Sichuan-Yunnan region from a number of focal mechanism solutions[J]. Acta Seismologica Sinica, 2006, 28(4): 337-346 (in Chinese with English abstract). Zhou Qingyun, Liu Zifeng, He Suge. Influence factors of velocity variation of airgun seismic waves and travel time changes related to earthquakes[J]. Earthquake, 2018, 38(3): 144-157 (in Chinese with English abstract). Zou Zhenxuan, Li Jinlong, Yu Tiehong. Determining average seismic wave velocity ratios (vP/vS) in Shanxi reservoir of Wenzhou region by using multi-station method[J]. Earthquake, 2006, 26(4): 133-137 (in Chinese with English abstract).