2. Office of the Western Yunnan Earthquake Prediction Study Area, China Earthquake Administration, Dali 671000, Yunnan, China
In the process of evolution of the earth's crust, the physical state of the medium will produce a series of changes, such as microfracturing, 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 normaldeclinerecoveryearthquake 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 reservoirinduced earthquakes, that is, the wave velocity ratio also showed a development process of a declinerecoveryearthquake (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 reservoirinduced earthquakes.
The above studies are all aimed at natural earthquakes or reservoirinduced 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 Pwave and Swave 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 Pwave and Swave 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 largevolume 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 YunlongEryuan 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 microdynamic change of physical properties of the underground medium. The paper also calculates the bvalue 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 PREPROCESSING 1.1 The Active Seismic Source Data of Binchuan Transmitting Seismic StationThe largevolume 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.
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 demeaning, detrending and antitilting to the airgun signals before picking up the phase data. Considering that the dominant frequency of airgun is concentrated at 25Hz (Chen Jia et al., 2016), the waveform is preprocessed by a 25Hz 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 Pwave 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 Pwave and Swave arrivaltimes 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 crosscorrelated 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 Pwave and Swave 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.
This paper uses the earthquake catalogue from January to June 2016 provided by the Yunnan Digital Seismic Network to calculate the bvalue of YunlongEryuan and its surrounding areas. In order to avoid the influence of data selection on the statistical results, the removal of the aftershock sequence, declustering of earthquake swarms and the analysis of completeness of the catalogue are crucial. This paper manually performed a second check based on the preliminary declustering 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, M_{C}, of the study area during the study period, which is M_{C}=1.0. The epicenter distribution of the earthquake with M_{L}≥1.0 in the study area from January to June 2016 is shown in Fig. 4.
The focal mechanism solutions used to calculate the consistency parameter are integrated by the relevant results obtained by the CAP method and the Pwave 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 M_{L}≥2.5 earthquakes in YunlongEryuan 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.
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 Pwave velocity is v_{P}, the travel time is t_{P}, Swave velocity is v_{S}, the travel time is t_{S}. 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, t_{S}=T_{S}T_{0}, t_{p}=T_{p}T_{0}; T_{p}, T_{S} are respectively the Pwave and Swave arrival time, T_{0} is the origin time of earthquake(in this paper, T_{0} 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, v_{P}/v_{S}, 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 T_{0}), as shown in the records of station EYA, the Pwave arrival time of this excitation is at 22:00:11. 983(namely T_{p}), the Swave arrival time is at 22:00:20.401(namely T_{S}), substituting T_{0}, T_{P} and T_{S} 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 timeseries variation characteristics. In order to reduce the effect of random variations on the calculation results, the SavitzkyGolay smoother integrated in Origin Software was used to smooth the calculation results of wave velocity ratio of each station in the timeseries variation analysis.
The SavitzkyGolay 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 SavitzkyGolay smoothing method is as follows:
Let a set of data in x(n) be x(i), i=M, …., 0…., M, by constructing a porder polynomial f_{i} 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_{k = 0}^p {{a_k}} {S_{k + r}} $  (6) 
Therefore, by selecting M (the number of singlesided points to be fitted) and p (polynomial order), F_{r} and S_{k+r} can be obtained, and substituting into equation (6), the coefficients a_{0}, a_{1}, …. a_{p} can be obtained. Thus, the polynomial f_{i} is determined, and the smoothed data value is obtained.
In this paper, the SavitzkyGolay smoother is used to smooth the wave velocity ratio of each station obtained by the singlestation 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.
GutenbergRicker'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 avalue represents the level of seismic activity in the statistical region, and the bvalue 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 bvalue 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 bvalue 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,
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 nonuniform stress fields in the focal mechanism solutions, Michael proposed a superposition stress field inversion method (SSI), which simulates the distribution of nonuniform stress fields by superimposing disturbances on a uniform stress field to obtain the spatial distribution and relative size φ=(S_{2}S_{3})/(S_{1}S_{3}) as well as the inversion variance of the three principal stress axes S_{1}, S_{2} and S_{3} (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 Moderatestrong EarthquakesThe wave velocity ratio before the M_{S}4.5 and M_{S}4.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 declinerecovery, 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 M_{S}4.5 and M_{S}4.0 earthquakes in Eryuan on February 8, 2016, the wave velocity ratio of stations 53036, EYA and YUL began to decrease quasisynchronously.
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 M_{S}5.0 and Eryuan M_{S}4.6 earthquakes on May 18, 2016. Since March 18, the wave velocity ratio of station 53036 experienced a change process of smallscale 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 declineearthquake; 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 declineturning 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 preearthquake 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 M_{S}5.0 and Eryuan M_{S}4.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 postearthquake stress adjustment.
Comparing the characteristics of the wave velocity ratios of the stations in different azimuths in the epicenter area before the M_{S}4.5 and M_{S}4.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 declinerecoveryrecovery speedup, but there were some differences in the time and range of the change. The DD mode mechanism of "expansionwater inflowing" can reasonably explain the development process of the declineriseearthquake of the preearthquake wave velocity, that is, the wave velocity decreases with the increase of fracture density in the dry fracture environment, and increases in a watersaturated 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, v_{P}/v_{S}, as well.
3.2 Analysis of the Characteristics of bvalue as a Function of TimeUsing the maximum likelihood method and the official catalogue provided by the Yunnan Digital Seismic Network with a threshold magnitude of M_{L}1.2 (>M_{C}), the daily sliding bvalue of accumulatively 15 days from January to June 2016 of the YunlongEryuan 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 bvalue 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 bvalue. According to the theory that the stress is inversely proportional to the bvalue (Scholz C. H., 1968; Wyss M. et al., 2000), it can be inferred that the decrease in the bvalue before the moderatestrong earthquake indicates the increase of stress, and the stress level decreases with the release of energy after an earthquake, so the bvalue began to rise after the earthquakes.
For the focal mechanism solutions of 57 M_{L}≥2.5 earthquakes in YunlongEryuan 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 nonuniform. 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.
In calculating the wavevelocity ratio using the singlestation method, the main factors affecting the calculation results are the pickup accuracy of arrival time of Pwave and Swave, 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 pickup accuracy of the arrival time of direct Pwave and Swave. In this paper, the waveform data is preprocessed before the arrival time is picked up by the waveform crosscorrelation method. Finally, with highquality 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^{ 3}s, the absolute arrival time difference error obtained is small, and the average standard deviation of the absolute arrivaltime difference of Pwave and Swave is 0.14 and 0.19, respectively. In addition, in analyzing the temporal variation characteristics of v_{P}/v_{S}, 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 v_{P}/v_{S} 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 DISCUSSIONBased on the excitation signals of the largecapacity airgun source in the Binchuan transmitting seismic station from January to June 2016, this paper acquired the absolute arrival time data of Pwave and Swave of the four stations near the epicenters of the Eryuan M_{S}4.5 and M_{S}4.0 earthquakes on February 8, 2016 and the Yunlong M_{S}5.0 and Eryuan M_{S}4.6 earthquakes on May 18, 2016, using the waveform crosscorrelation method. The wavevelocity ratios of the four stations are calculated by the singlestation method, and at the same time, the bvalue, the focal mechanism consistency parameter, and variance of YunlongEryuan 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 declinerecoveryfast rise, and the time and range of the change were different. The quasisynchronization of rapid rise occurring within 312 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 M_{S}4.5 and M_{S}4.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 bvalue and the focal mechanism consistency parameter, variance, commonly used to characterize the stress level were calculated. The bvalue 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, bvalue 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.
ACKNOWLEDGEMENTThe author would like to express their heartfelt thanks to all the institutions and individuals who have contributed to the largevolume airgun experiment in Binchuan.
Aki K. Maximum likelihood estimate of b in the formula logN=abM and its confidence limits[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 1965, 43(2): 237239. 
Cai Tianjing, Tang Han. Review on the leastsquaresfit principle of SavitzkyGolay smoothing filter[J]. Digital Communication, 2011, 38(1): 6368, 82 (in Chinese with English abstract). 
Chen Jia, Li Xiaobin, Yang Jun, Ye Beng. Large volume airgun source spectrum characteristics of Binchuan, Yunnan[J]. Earthquake Research in China, 2016, 32(2): 216221 (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 M_{S}5.0 earthquake using airgun signals[J]. Journal of Seismological Research, 2017, 40(4): 550556, 677 (in Chinese with English abstract). 
Chen Yang, Lv Yuejun, Xie Zhuojuan, Pan Long. Review of the study of seismicity parameter bvalue[J]. Bulletin of the Institute of Crustal Dynamics, 2013, 3847 (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): 142159 (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): 17391751. DOI:10.1007/s1143001690966 
Feng Deyi, Tan Aina, Wang Kefen. Velocity anomalies of seismic waves from near earthquakes and earthquake prediction[J]. Acta Geophysica Sinica, 1974, 17(2): 8498 (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): 292309 (in Chinese with English abstract). 
Feng Deyi. Velocity Anomaly of Seismic Waves[M]. Beijing: Seismological Press, 1981. 215, 140149, 236246 (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): 93059320. DOI:10.1029/JB089iB11p09305 
Hardebeck J.L., Michael A.J. Damped regionalscale 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. Spatialtemporal characteristics of the focal mechanism consistency parameter in Tianshan (within Chinese territory) seismic zone[J]. Seismology and Geology, 2015, 37(3): 792803 (in Chinese with English abstract). 
Li Mingxiao, Liu Jie. Study on velocity ratio (v_{P}/v_{S}) anomaly of earthquake sequences in Yunnan region[J]. Earthquake, 2006, 26(1): 2634 (in Chinese with English abstract). 
Li Yan'e, Wang Linying, Zheng Xuyao. Restudy of the variation of v_{P}/v_{S} before and after the Wenchuan earthquake[J]. Acta Seismologica Sinica, 2014, 36(3): 425432 (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): 2530 (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): 464470 (in Chinese with English abstract). 
Liu Xuzhou. Study on Seismic TravelTime 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 M_{S}5.0 earthquake in 2018[J]. Journal of Seismological Research, 2018, 41(4): 494502 (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): 53835402. 
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): 213217 (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): 1151711526. 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): 357368. 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): 14411444. 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): 63036319. DOI:10.1029/91JB00195 
Scholz C.H. The frequencymagnitude relation of microfracturing in rock and its relation to earthquakes[J]. Bulletin of the Seismological Society of America, 1968, 58(1): 399415. 
Silver P.G., Daley T.M., Niu F.L., Majer E. L. Active source monitoring of crosswell seismic travel time for stressinduced changes[J]. Bulletin of the Seismological Society of America, 2007, 97(1B): 281293. 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 bvalue before major earthquakes in the SichuanYunnan region[J]. Seismology and Geology, 2012, 34(3): 531543 (in Chinese with English abstract). 
Wang Linying, Guo Yongxia, Liu Fang, Jiang Changsheng. Temporal v_{P}/v_{S} variation characteristics in different zones of China's Capital Circle area before and after Wen'an earthquake[J]. Acta Seismologica Sinica, 2008, 21(3): 243257 (in Chinese with English abstract). DOI:10.1007/s115890080243z 
Wang Linying, Li Yan'e, Zheng Xuyao, Wang Shengwen. Temporal variation of v_{P}/v_{S} at single seismic station before the 2013 Lushan M_{S}7.0 earthquake[J]. Acta Seismologica Sinica, 2014, 36(1): 4258 (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): 223233 (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 airgun in Xinjiang before two earthquakes with M_{S}5.0[J]. Earthquake Research in China, 2016, 32(2): 270281 (in Chinese with English abstract). 
Weng Zhaoqiang, Liang Xiangjun, Wu Haoyu, Liu Linfei, Li Li. The variation of v_{P}/v_{S} at single seismic station of Shanxi region studied by Wadati method of single station and multiearthquake[J]. Earthquake Research in China, 2015, 31(1): 89100 (in Chinese with English abstract). 
Wiemer S. A software package to analyze seismicity:ZMAP[J]. Seismological Research Letters, 2001, 72(2): 373382. 
Wiemer S., Gerstenberger M., Hauksson E. Properties of the aftershock sequence of the 1999 M_{W}7.1 Hector Mine earthquake:Implications for aftershock hazard[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 12271240. DOI:10.1785/0120000914 
Wyss M., Schorlemmer D., Wiemer S. Mapping asperities by minima of local recurrence time:San JacintoElsinore fault zones[J]. Journal of Geophysical Research:Solid Earth, 2000, 105(B4): 78297844. 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): 715 (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): 476486 (in Chinese with English abstract). 
Yi Guixi, Wen Xueze, Su Youjin. Study on the potential strongearthquake risk for the eastern boundary of the SichuanYunnan active faultedblock, China[J]. Chinese Journal of Geophysics, 2008, 51(6): 17191725 (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): 282293 (in Chinese with English abstract). 
Yi Guixi, Wen Xueze, Xin Hua, Qiao Huizhen, Wang Siwei, Gong Yue. Stress state and majorearthquake risk on the southern segment of the Longmen Shan fault zone[J]. Chinese Journal of Geophysics, 2013, 56(4): 11121120 (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. Traveltime variations before and after two major earthquakes derived from activesource seismic data[J]. Chinese Journal of Geophysics, 2017, 60(10): 38153822 (in Chinese with English abstract). 
Zhang Zhiwei, Zhou Longquan, Long Feng, Ruan Xiang. Spatial and temporal characteristic of stress field for Wenchuan M_{S}8.0 and Lushan M_{S}7.0 earthquake sequence[J]. Seismology and Geology, 2015, 37(3): 804817 (in Chinese with English abstract). 
Zhong Jimao, Cheng Wanzheng. Determination of directions of the mean stress field in SichuanYunnan region from a number of focal mechanism solutions[J]. Acta Seismologica Sinica, 2006, 28(4): 337346 (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): 144157 (in Chinese with English abstract). 
Zou Zhenxuan, Li Jinlong, Yu Tiehong. Determining average seismic wave velocity ratios (v_{P}/v_{S}) in Shanxi reservoir of Wenzhou region by using multistation method[J]. Earthquake, 2006, 26(4): 133137 (in Chinese with English abstract). 