Seismic wave is the most effective way for us to understand the structure of the earth and the state of the earths medium. Most of our understanding of the earths interior comes from seismic wave. According to different research objectives and contents, researches using seismic waves can be roughly divided into underground medium structure and medium properties change with time. In medium structure research, signaltonoise ratio (SNR) of seismic signals is improved by means of encrypted observation array and longtime superposition, so as to improve the spatial resolution and accuracy of the structure research results (Yao Huajian et al., 2010; Bai Zhiming et al., 2016). The study of medium properties change with time focuses on the changing trend of the properties of underground media with time and the relationship between them and earthquakes (Chen Yong et al., 2016). Besides high signaltonoise, high repeatability of seismic signals at different times is also needed in medium properties change reasearch.
The stress in the crustal medium accumulates slowly before the earthquake, releases sharply during the earthquake and rises slowly after the earthquake. Under different stress states, the pore in the crustal medium will close or open, and the wave velocity and anisotropic of fast wave direction of the medium will change accordingly. Some scholars use data from passive sources such as seismic waves from natural earthquakes and noise to research medium properties change (Poupinet G. et al., 1984; Brenguier F. et al., 2008; Cheng Xin et al., 2010; Hillers G. et al., 2015). However, passive sources exist some natural defects, such as the uncertainty of the spatial location of natural earthquakes and the low temporal resolution of noise data. In recent years, a type of new active source called airgun used in the ocean has been introduced to form an efficient and highquality artificial source in terrestrial reservoirs by Chinese scholars (Chen Yong et al., 2017; Wang Baoshan et al., 2016, 2018). Compared with previous artificial sources such as hammering and blasting, the airgun source has the advantages of high repeatability, long propagation distance, and no damage to the excitation environment, and can greatly improve the signaltonoise ratio and propagation distance by superposition. These characteristics help us to improve the accuracy and time resolution in the study of crustal media.
Since 2012, some earthquake signal transmitting seismic stations using airgun have been built in China, such as Binchuan, Yunnan, Hutubi, Xinjiang, Qilian Mountain, Gansu, etc. Using the data accumulated by these airgun source, scholars have done many research works(Wang Bin et al., 2016; Zhang Yuansheng et al., 2016; Su Jinbo et al., 2016). She Yuyang et al. (2019) studied the characteristics of seismic wave propagation of Binchuan airgun signal. Su Jinbo et al. (2019) studied the seasonal variation of Hutubi airgun signal. Yang Wei et al. (2016) and Chen Huifang et al. (2016) studied the effects of excitation environment on airgun source. Wei Congxin et al. (2018) studied the optimal depth of airguns. Tian Xiaofeng et al. (2018) studied the crustal velocity structure of the Yangtze River in Anhui using Pwave travel time tomography of airgun signals. She Yuyang et al. (2018) generated the shallow crustal velocity structure in the Anhui section of the Yangtze River using surface wave of airgun signals. These studies have deepened our understanding of the airgun source and facilitated further research related to airguns (Wang Weitao et al., 2017). Some useful data processing methods of airgun signal have been proposed by some scholars also (Tan Junqing et al., 2019; Li Jun et al., 2019; Wang Weijun et al., 2019).
In the study of changes of crustal medium by airgun sourceand other active source, the main research object is medium wave velocity, including Pwave velocity, Swave velocity, coda wave velocity, etc. For example, Liu Zifeng et al. (2019) used the data of airgun source to study the Wave Velocity Ratio of airgun signals before and after two major earthquakes. Yang Jianwen et al. (2019) used airgun source signals to study regional wave velocity changes before and after the Yunlong M_{S}5.0 and Yangbi M_{S}5.1 Earthquakes. Yang Wei et al. (2018) studied the influence of groundwater level on the study of the travel time variation of airgun signals. These studies show that the airgun source can be used to monitor the stress changes in crustal media. Based on the data from November to December, 2015 at BESTSS, this paper attempts to study the influence of crustal medium stress on the direction of airgun signal propagation, expecting to obtain another physical quantity to monitor the stress changes in crustal media besides the variation of velocity.
1 DATAThe Binchuan Earthquake Signal Transmitting Seismic Station (BESTSS) in Yunnan Province was built in 2012, and a rectangular array of four airguns was suspended 10m below the water surface of Dayindian reservoir. Since the completion of the construction, except for the period of low water level of the reservoir from May to August in each year, airgun array shotted every other day. Studies have shown that when the airgun source is stopped for a long time and then reused, the airgun signals repetition rate is low, and changes in reservoir water level, excitation pressure and other factors will also affect the signal repetition rate (Zhou Qingyun et al., 2018). Based on the above reasons, we selected waveform of 943 shots intensively excited from November 18 to December 13, 2015. During the excitation period, the water level of reservoir changed little, the excitation pressure was fixed, and the repeatability of airgun signals was good. During the excitation period, 012 oclock every hour per shot, 1224 oclock every 15 minutes per shot, and from November 29 and December 1 due to airgun maintenance, no excitation.
During the experiment period, there were 48 mobile seismic stations, with the farthest distance of 151km and the nearest 0.8km (Fig. 1). The sampling rate of all stations was 100Hz. We study the change of azimuth angle with time after waveform propagation, so we need a higher signaltonoise ratio of each waveform. Therefore, we use the data of stations with less noise and clear signal of each shot.
We analyzed the repeatability of the signals excited by the airgun source during the experiment using CKT1 (epicenter distance ~40m) on the shore of the reservoir. The correlation coefficients greater than 0.995 accounted for 94.3% and that greater than 0.999 accounted for 40.4%. It shows that the signals generated by the airgun source during the experiment have excellent repeatability (Fig. 2). The maximum amplitude of the airgun signal during the test is also analyzed. It is 88.2% in the range of average amplitude (±1%) and 99.3% in the range of (±2%). It also shows that the excitation condition of the airgun has not changed during the test (Fig. 2).
Stockwell R. G. et al., (1996) proposed the Stransform method, which is a timefrequency analysis method. This method combines the advantages of shorttime Fourier transform and wavelet transform. Based on Stransform, scholars have proposed different filtering methods. Zheng Chenglong et al. (2015) proposed Stransform template filtering based on the high repeatability of airgun signal.
Linear superposition has minimal waveform distortion (Wu Anxu et al., 2016), so we linearly superimpose the 943 signals when creating the filter template. After stacking, each station signal has a high signaltonoise ratio. The timefrequency spectrum S(τ, f) of the superimposed signal h(t) is obtained by Stransform, and the timefrequency spectrum is normalized. Set a threshold a, when the value in the spectrum S(τ, f) is smaller than a, the filter template D(τ, f)=0, otherwise D(τ, f)=1. The result of seismic signal filtering is obtained by Stransform, multiplication with filter template and Sinverse transformation into time domain.
2.2 Calculation of Back Azimuth of Airgun SignalThe Pwave of airgun signal propagates to seismic station through crustal medium. Its particle trajectory should be roughly in the plane including source and receiving station. If we rotated the horizontal component recorded by seismograph to the radial and tangential direction, the projection of Pwave in the horizontal plane mainly concentrates on the radial component, and the energy in the tangential component is very small. To calculate the back azimuth angle, the covariance matrix of two horizontal components should be first calculated (Niu Fengling et al., 2011).
$ c_{i j}=\int_{0}^{T} u_{i}(t) u_{j}(t) \mathrm{d} t \quad i, j=1, 2 $  (1) 
Where T is the Pwave window length, and u_{i}(t) and u_{j}(t) are the horizontal components of the seismic record. In the absence of noise, the covariance matrix c has a nonzero eigenvalue, and the direction of the eigenvector corresponding to the eigenvalue is the back azimuth angle:
$ \theta=\cot ^{1} \frac{c_{22}c_{11}+\sqrt{\left(c_{11}c_{22}\right)^{2}+4 c_{12}^{2}}}{2 c_{12}} $  (2) 
When there is noise, c will have two nonzero eigenvalues, and the ratio of the two nonzero eigenvalues reflects the linearity of the particle motion. The Stransformation template filter can significantly suppress the noise in the airgun signal, and the signaltonoise ratio can be used to estimate the error in the azimuth calculation (Drew J. et al., 2008).
2.3 Data Processing FlowFirst, preprocessing. Waveforms taken from continuous waveform recordings typically contain nonzero constants, linear trends, and noise, so proper preprocessing is required. The preprocessing mainly consists of three steps: detrending, removing the linear component of the waveform record; deaveraging, removing the nonzero constant of the waveform; and bandpass filtering to initially reduce some noise. Considering the main frequency range of the airgun source, the bandpass filtering range in the preprocessing step is 110Hz, which removes some common noises and ensures that the airgun signal is not filtered.
Then, filtering using Stransform template filter. The airgun signal recorded by stations located 10km away from the epicenter is relatively low in signaltonoise ratio and cannot be used directly. Stransform template filtering can suppress most noise under the premise of retaining effective airgun signal, and effectively improve the signaltonoise ratio. The template filtering mainly includes several steps: linearly superimposing the airgun signals of the station, performing S transformation and normalization; selecting appropriate thresholds to create a filtering template; performing S transformation on each signal and multiplying it by the filtering template. Then inversely Stransform to the time domain. The filter template threshold is 0.15. It can be seen from the filtering results that the Stransform template filtering can greatly suppress the noise while effectively retaining the target signal (Fig. 3).
Finally, selecting appropriate time window of the horizontal component of each stationt to calculate the azimuth of the airgun signal. In the selection of the time window, the starting time is the time when the Pwave arrives at the station (Pwave velocity is about 5.55km/s in research aera). The choice of termination time should consider the influence of the duration of the focal time function on the nearsource stations and the arrival time of the Swave at the farreaching stations. Therefore, we use the termination time as follows:
$ {t_{{\rm{end }}}} = {t_{{\rm{porrive }}}} + {\rm{dist}}/14 + 3.25 $  (3) 
Where t_{p_arrive} is the time when the Pwave arrives at the station, and dist is the epicenter distance of the station.
3 DATA PROCESSING RESULTSWhen no earthquake occurs, the stress in the crustal medium does not change much, and the change of the back azimuth angle is small. In fact, due to the influence of signaltonoise ratio and the condition of stations foundation, the back azimuth angle changes of 943 shots are relatively large in some stations. We consider that the reliability of the back azimuth angle calculated from the data of these stations is insufficient. The stations with an average angle deviation of less than 5 degrees are considered to be of good quality. A total of 22 stations meet this condition. The spatial distribution of the stations is shown in Fig. 1, and the calculation results of the back azimuth are shown in Fig. 4 and Fig. 5. According to the result, we can obtain two preliminary conclusions:
Firstly, the results in Fig. 3 show a diurnal variation characteristic. By enlarging result pictures of some stations in Fig. 4, it can be seen that the change of back azimuth angle has clear diurnal variation characteristics (Fig. 6). The regular change in 10 consecutive days indicates that the diurnal variation characteristics are not caused by the accidental factors such as the position of airgun and the excitation conditions. It is presumed that the diurnal variation may be caused by the change of tidal stress. Sensitivity of physical parameters response to stress changes is an important index to determine whether a physical parameter can be used to detect stress changes in crustal media. For example, the wave velocity variation and travel time variation of the crustal medium are highly sensitive to stress changes like tidal stress change. Both of them change significantly, and the trend is highly correlated with tidal stress change (Niu Fengling et al., 2008; Wang Baoshan et al., 2008), so these two physical parameters can be used to monitor stress changes in the crustal medium. In this study, the azimuth of seismic ray is sensitive to the change of tidal stress, so it may be used as a physical parameter to monitor the change of stress in crustal medium.
Secondly, the azimuth angle of 12 stations in Fig. 4 varies greatly from 16:00 to 24:00 on December 4, ranging from 6° to 10°. Stations with larger azimuth change, smaller azimuth change and no obvious change have no clear distribution patterns in space (Fig. 1). Two stations not far away from each other even show opposite trends, such as 53273 and EY01, 53274 and EY02.
4 DISCUSSION 4.1 Error AnalysisIn the absence of noise, the azimuth calculated by equation (2) is the direction of Pwave vibration; when there is noise, the calculated azimuth may not coincide with the direction of Pwave vibration (Niu Fengling et al., 2011). In order to evaluate the influence of noise on the azimuth calculation, we add real noise received by the station with different intensity (gain) to the superposition of airgun signal received in station EY02, and calculate the azimuth under different SNR. The angle is compared to the azimuth calculation without adding noise (Fig. 7). The station EY02 received a total of 823 airgun signals, and we selected 2 segments of noise before each signal, using a total of 1646 segments of noise data to calculate the impact of larger and smaller noise. We performed 40 kinds of gain level on each noise segment, from 10^{0.1} to 10^{4}. We found that the higher the signaltonoise ratio of the airgun signal, the smaller the error. When the SNR>10, the error σ < 1°, when the SNR>80, the error σ < 0.1°, which is similar to the research results of other scholars (Drew J. et al., 2008; Moriya H. et al., 1994; Grechka V. et al., 2011).
After Stransform template filtering, the amplitude of the noise portion is set to 0, so the traditional method cannot be used to calculate the signaltonoise ratio of filtered airgun signal. The method used in this paper is: select 10 segments of continuous noise records with 20s length randomly in the first 100 seconds of Pwave arrival, filter the noise, and use the filtered data to evaluate the noise level; after Pwave arrival of the filtered airgun signal, the data within 20 seconds are used to evaluate the signal level. The formula for calculating SNR is as follows:
$ {\mathop{\rm SNR}\nolimits} = \frac{{\sum\limits_{j = 1}^{2000} {{S_j}} \times {S_j}}}{{\frac{1}{{10}}\sum\limits_{n = 1}^{10} {\sum\limits_{i = 1}^{2000} {{N_{ni}}} } \times {N_{ni}}}} $  (4) 
Where S_{j} is the value of the jth sampling point of the airgun signal, and N_{ni} is the value of the ith sampling point of the nth segment of noise data.
We calculated the signaltonoise ratio of each airgun signal for the 22 stations used in this paper, and the results are shown in Fig. 8. It can be seen from the result graphs that the signaltonoise ratios of stations within 10km of the airgun source are greater than 80. The error of the azimuth calculated by these stations is less than 0.1°; the signaltonoise ratios of stations above 10km are mostly greater than 10, and the azimuth error is less than 1°. The diurnal variation of azimuth angle is about ±0.5°±1°, the amplitude of the azimuth change of all stations on December 4th is ≥2°, both of them are greater than the influence of noise, therefore, the calculated results of azimuth change are reliable.
There are pores in the rock. When the intensity or direction of the stress in the rock changes, the openclose degree and direction of the rock pores will change (Crampin S., 1987), in which the change of pores openclose degree leads to the change of wave velocity (Wang Baoshan et al., 2008; Niu Fenglin et al., 2008), and changes of direction of pores will cause changes of anisotropic of fast wave direction of the rock (Montagner J. P. et al., 1986; Durand S. et al, 2011).
The rock types in research area mainly consist of clastic rocks, limestone, mudstone and sandstone. The sensitivity of these rocks to stress varies by more than 13 times (Birch F., 1960, 1961; Christensen N. I., 1965; Freund D., 1992). When the underground stress conditions change, the velocity change of different rocks varies. If the Pwave obliquely enters the interface of different rocks, the wave propagation direction will change. According to Snells law, the change of the propagation direction is given by:
$ \Delta \theta_{1}=\arcsin \left(\frac{v_{1}\left(1+\Delta P * \alpha_{1}\right)}{v_{2}\left(1+\Delta P * \alpha_{2}\right)}\right)\arcsin \left(\frac{v_{1}}{v_{2}} * \sin \theta_{2}\right) $  (5) 
Where v_{1} and v_{2} are the Pwave velocities of the media on two sides of the interface, α_{1} and α_{2} are the stress sensitivities, ΔP is the maximum tidal stress change, and θ_{2} is the incident angle. Take v_{1}=3.0km/s, v_{2}=4.5km/s, α_{1}=1.2×10^{－6}Pa^{－1}, α_{2}=4×10^{－7}Pa^{－1}, ΔP=5kPa. When the incident angle θ_{2} is between 15°40°, the Δθ_{1} value is 0.04°0.10°, which is about 10 times smaller than the previous results. Therefore, it is considered that the difference in wave velocity variation of different rocks is not the cause of Pwave azimuth change.
A study carried by Durand S. et al., (2011) showed that when the velocity of crustal medium with 3% anisotropy changes by 1‰, the azimuth change after propagation is about 2.5°. The daily wave velocity change calculated by the data of the Binchuan airgun signal monitoring network is about 1‰, and the azimuth change calculated in this paper is about 1°2°, and the consistency between the two is good.
The results in Fig. 8 shows that the signaltonoise ratio of airgun signal has clear diurnal variation characteristics. Because the energy variation range of airgun signal is very small during this period, the diurnal variation characteristic response that noise has diurnal characteristics: human activities are frequent during the day, so the noise intensity is high and the signaltonoise ratio is low; human activities are weakened at night, reduced noise intensity, and increased signaltonoise ratio. We believe that the diurnal variation of azimuth does not come from the daily variation of noise. The reasons are: ① the direction of motion of the particle caused by noise is random, and its influence on the calculation result of azimuth is increasing of error; ② The Stransform template filtering method can suppress the noise very well; ③ if the azimuth variation feature is derived from noise, the farther the epicenter distance and the lower the signaltonoise ratio are, the azimuth calculation result should have more obvious diurnal characteristics—this is contrary to our calculated results. Therefore, we believe that the azimuthal variation feature does not come from the daily variation of noise.
According to the comprehensive analysis, the diurnal variation of the azimuth angle of the Pwave of the airgun signal is mainly caused by the change of the anisotropic of fast wave direction caused by the stress change.
4.3 Azimuthal Abrupt Change Analysis on December 4The details of the azimuth change on December 4 cannot be seen clearly in Fig. 4. Therefore, the azimuth change map of 12 stations in Fig. 4 and the correlation coefficient map in Fig. 2 are enlarged and placed in Fig. 9. It can be seen from Fig. 9 that at 17:30 on December 4, the correlation coefficient of the airgun signal began to decrease, and the azimuth of the station began to change; after 0:00 on December 5, the correlation coefficient of the airgun signal suddenly returned to normal(corr > 0.999), and the azimuth of the airgun signal of each station also suddenly recovered. There was no reservoir water level, temperature, etc. during this period, and there was no rain or earthquake bigger than magnitude 1.0. Therefore, the reason for the sudden change in azimuth on December 4 was the sudden change of the airgun source.
One of the reasons why the airgun source is better than natural earthquakes and noise is that it can be frequently excited and has a high time resolution. For stations far away from epicentres, the amplitude of singleshot signal is very small, and it is often necessary to superimpose multiple times to obtain a reliable signal, but the superposition will reduce the time resolution. In this paper, the Stransform template filtering method is used to filter the singleshot signal, and then the azimuth is calculated. This method also has better results in stations with a slightly longer epicenter. We selected two representative Station 53251 and 53252 which far away from airgun. Their epicenter distances were 87km and 151km respectively. The calculation results of the azimuth changes of the two stations are shown in Fig. 10. The direct calculation results are more discrete and have no obvious pattern. In order to improve the stability of the calculated results, we have adopted two processing methods: every filtered airgun signal is superimposed with its adjacent 4 shots and the azimuth is calculated; the result of the single shot calculation is smoothed every 5 adjacent points. It can be seen from Fig. 10 that after superposition or smoothing, the back azimuth of the station with a farther distance from the epicenter has better stability and daily variation.
The above calculation results are obtained when the threshold is 0.15. In order to study the influence of the filtering template threshold selection on the calculation result of the back azimuth change, we changed the threshold from 0.15 to 0.05 and 0.30, and recalculated the change of the back azimuth. The results show that the calculation results of the filtering templates of three different thresholds are highly similar. We calculated the correlation coefficient of filtered waveform of the same station between template threshold 0.05 and template threshold 0.15, template threshold 0.30 and template threshold 0.15, and the correlation coefficient calculated by all stations is plotted in Fig. 11 according to the epicenter distance. It can be seen that stations with an epicenter distance of 151km have a correlation coefficient of more than 0.97, and the correlation coefficient of stations within 25km exceeds 0.99. The reason for this phenomenon is that the superimposition can suppress the noise and amplification airgun signal, and the superposition of 943 waveforms has already suppressed the noise to below 0.05 and amplified the airgun signal to 0.30 or more (in normalized spectrum). The above analysis shows that the selection of the threshold of the filter template in this study has little effect on the results, even the threshold value is changed from 0.05 to 0.30, the same result can be obtained.
The calculated azimuth of the airgun signal is different from the real azimuth of the station. There may be three causes:
First, the selection of the time window calculated the azimuth. Ideally, a window containing only Pwave phase should be the best selected. Since the source time function of the Binchuan active source signal lasts longer than 3s, the Pwave and Swave phases will separate only when the epicenter distance >10km. This may be the main reason for the deviation between the calculated azimuth and the true azimuth of the station within 10km of the epicenter. The second cause is the impact of complex shallow media. The airgun source has the characteristics of surface excited surface receiving and high frequency. It is sensitive to the difference of wave velocity in different media in the shallow crust. There are Quaternary loose sediments with great slow wave velocity (such as Binchuan Basin), soft lowspeed rock (such as clastic rock Q_{p}, mudstone J_{3}) and hard highspeed rock (such as limestone D_{1}, basalt P_{β}) medium in research area, so airgun signals do not propagate along a straight line, but bypass some basins with lower wave speeds or propagated along the faster media curves. This may be the main reason for the azimuth deviation of longdistance stations. Third, the stations in the active source network are all mobile stations, so the north of the instrument may have a deviation of several degrees to more than a dozen degrees from the geographic north.
Our purpose in this paper is the change trend of azimuth angle with time after airgun signal propagation. The difference between real azimuth angle and calculated azimuth angle has little influence on the research results, so no further analysis is made here.
4.7 The Influence of Stransform Template Filtering on Azimuth CalculationFiltering, superposition and other methods may cause waveform distortion (Wu Anxu et al., 2016). The filtering method used in this paper may also cause waveform distortion and affect the azimuth calculation results. We selected 53273 and EY01 stations with nearer epicenter distance, and calculated the azimuth variation and the correlation coefficient between each signal and superimposed waveform by Stransform template filtering and nonStransform template filtering. The result is shown in Fig. 12. It can be seen that the filtered waveform has a higher correlation coefficient; the filtering has a similar tendency to the azimuth calculated for the filtering. The above results show that the Stransform template filtering method can improve the waveform quality and stability, and has little influence on the azimuth variation analysis.
(1) The back azimuth of the airgun signal propagation is sensitive to stress changes. We calculated the propagation direction of the airgun signal received by the station and found that some stations have a diurnal variation characteristic. This characteristic is due to the change of tidal stress. The change of tidal stress is about 5kPa in one day. The stress change caused by moderatestrong earthquake is generally much larger than this value. Therefore, the azimuth angle of the airgun signal can be used to monitor the change of the stress state of the crustal medium.
(2) At the afternoon of December 4, the azimuth of airgun signal at more than half of the stations changed abruptly. After eliminating the possible influencing factors such as water level, temperature and precipitation of the reservoir, according to the consistency between the correlation coefficient of airgun signal and the azimuth change, it was considered that the abrupt change of azimuth resulted from the change of airgun source.
(3) When the receiving station is far away, the propagation azimuth stability calculated by the Stransform template filtered singleshot signal is poor. The 5shot superposition or 5point smoothing method can effectively improve the stability of the calculated results.
Bai Zhiming, Wu Qingju, Xu Tao, Wang XiaoBai Zhiming, Wu Qingju, Xu Tao, Wang Xiao. Basic features of the crustal structure in the lower Yangtze and its neighboring area in the Chinese Mainland:review of deep seismic sounding research[J]. Earthquake Research in China, 2016, 30(3): 298315. 
Birch F.Birch F. The velocity of compressional waves in rocks to 10 kilobars:1[J]. Journal of Geophysical Research, 1960, 65(4): 10831102. DOI:10.1029/JZ065i004p01083 
Birch F.Birch F. The velocity of compressional waves in rocks to 10 kilobars:2[J]. Journal of Geophysical Research, 1961, 66(7): 21992224. DOI:10.1029/JZ066i007p02199 
Brenguier F., Campillo M., Hadziioannou C., Shapiro N.M., Larose E.Brenguier F., Campillo M., Hadziioannou C., Shapiro N.M., Larose E. Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations[J]. Science, 2008, 321(5895): 14781481. DOI:10.1126/science.1160943 
Chen Huifang, Lin Binhua, Jin Xing, Wu Lihua, Cai HuitengChen Huifang, Lin Binhua, Jin Xing, Wu Lihua, Cai Huiteng. An experimental study on optimization of largevolume airgun source excitation conditions in a reservoir[J]. Earthquake Research in China, 2016, 30(3): 355363. 
Chen Yong, Wang BaoshanChen Yong, Wang Baoshan. Preface to the special issue on advances in seismic exploration and monitoring with active sources in China[J]. Earthquake Research in China, 2016, 30(3): 281283. 
Chen Yong, Wang Baoshan, Yao HuajianChen Yong, Wang Baoshan, Yao Huajian. Seismic airgun exploration of continental crust structures[J]. Science China Earth Sciences, 2017, 60(10): 17391751. DOI:10.1007/s1143001690966 
Cheng Xin, Niu Fenglin, Wang BaoshanCheng Xin, Niu Fenglin, Wang Baoshan. Coseismic velocity change in the rupture zone of the 2008 M_{W}7.9 Wenchuan earthquake observed from ambient seismic noise[J]. Bulletin of the Seismological Society of America, 2010, 100(5B): 25392550. DOI:10.1785/0120090329 
Christensen N.I.Christensen N.I. Compressional wave velocities in metamorphic rocks at pressures to 10 kilobars[J]. Journal of Geophysical Research, 1965, 70(24): 61476164. DOI:10.1029/JZ070i024p06147 
Crampin S.Crampin S. Geological and industrial implications of extensivedilatancy anisotropy[J]. Nature, 1987, 328(6130): 491496. DOI:10.1038/328491a0 
Drew J., White R., Wolfe J. Microseismic event azimuth estimation: establishing a relationship between hodogram linearity and uncertainty in event azimuth. In: Proceedings of 2008 SEG Annual Meeting. Las Vegas, Nevada: SEG, 2008: 14461450.

Durand S., Montagner J.P., Roux P., Brenguier F., Nadeau R.M., Ricard Y.Durand S., Montagner J.P., Roux P., Brenguier F., Nadeau R.M., Ricard Y. Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake[J]. Geophysical Research Letters, 2011, 38(13): L13303. DOI:10.1029/2011GL047875 
Freund D.Freund D. Ultrasonic compressional and shear velocities in dry clastic rocks as a function of porosity, clay content, and confining pressure[J]. Geophysical Journal International, 1992, 108(1): 125135. DOI:10.1111/j.1365246X.1992.tb00843.x 
Grechka V., Singh P., Das I.Grechka V., Singh P., Das I. Estimation of effective anisotropy simultaneously with locations of microseismic events[J]. Geophysics, 2011, 76(6): WC143WC155. DOI:10.1190/geo20100409.1 
Hillers G., BenZion Y., Campillo M., Zigone D.Hillers G., BenZion Y., Campillo M., Zigone D. Seasonal variations of seismic velocities in the San Jacinto fault area observed with ambient seismic noise[J]. Geophysical Journal International, 2015, 202(2): 920932. DOI:10.1093/gji/ggv151 
Li Jun, You Xiuzhen, Hu ShufangLi Jun, You Xiuzhen, Hu Shufang. Probability data screening method in airgun signal processing application[J]. Earthquake Research in China, 2019, 33(2): 235247. 
Liu Zifeng, Ye Beng, Chen Jia, Zhou Qingyun, Wang Guangming, Peng Guanling, Li YongliLiu Zifeng, Ye Beng, Chen Jia, Zhou Qingyun, Wang Guangming, Peng Guanling, Li Yongli. Study of the change in wave velocity ratio before and after two strong earthquakes using airgun source data[J]. Earthquake Research in China, 2019, 33(2): 305319. 
Montagner J.P., Nataf H.C.Montagner J.P., Nataf H.C. A simple method for inverting the azimuthal anisotropy of surface waves[J]. Journal of Geophysical Research, 1986, 91(B1): 511520. DOI:10.1029/JB091iB01p00511 
Moriya H., Nagano K., Niitsuma H.Moriya H., Nagano K., Niitsuma H. Precise source location of AE doublets by spectral matrix analysis of triaxial hodogram[J]. Geophysics, 1994, 59(1): 3645. DOI:10.1190/1.1443532 
Niu Fenglin, Silver P.G., Daley T.M., Cheng Xin, Majer E.L.Niu Fenglin, Silver P.G., Daley T.M., Cheng Xin, Majer E.L. Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site[J]. Nature, 2008, 454(7201): 204208. DOI:10.1038/nature07111 
Niu Fenglin, Li JuanNiu Fenglin, Li Juan. Component azimuths of the CEArray stations estimated from Pwave particle motion[J]. Earthquake Science, 2011, 24(1): 313. 
Obermann A., Froment B., Campillo M., Larose E., Planès T., Valette B., Chen J. H., Liu Q. Y. S.Obermann A., Froment B., Campillo M., Larose E., Planès T., Valette B., Chen J. H., Liu Q. Y. S. eismic noise correlations to image structural and mechanical changes associated with the M_{W}7.9 2008 Wenchuan earthquake[J]. Journal of Geophysical Research, 2014, 119(4): 31553168. DOI:10.1002/2013JB010932 
Poupinet G., Ellsworth W.L., Frechet J.Poupinet G., Ellsworth W.L., Frechet J. Monitoring velocity variations in the crust using earthquake doublets:An application to the Calaveras Fault, California[J]. Journal of Geophysical Research, 1984, 89(B7): 57195731. DOI:10.1029/JB089iB07p05719 
She Yuyang, Yao Huajian, Zhai Qiushi, Wang Fuyun, Tian XiaofengShe Yuyang, Yao Huajian, Zhai Qiushi, Wang Fuyun, Tian Xiaofeng. Shallow crustal structure of the middlelower Yangtze river region in eastern China from surfacewave tomography of a large volume airgunshot experiment[J]. Seismological Research Letters, 2018, 89(3): 10031013. DOI:10.1785/0220170232 
She Yuyang, Yao Huajian, Wang Weitao, Liu BinShe Yuyang, Yao Huajian, Wang Weitao, Liu Bin. Characteristics of seismic wave propagation in the Binchuan region of Yunnan using a dense seismic array and large volume airgun shots[J]. Earthquake Research in China, 2019, 33(2): 174185. 
Stockwell R.G., Mansinha L., Lowe R.P.Stockwell R.G., Mansinha L., Lowe R.P. Localization of the complex spectrum:the S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 9981001. DOI:10.1109/78.492555 
Su Jinbo, Wang Qiong, Wang Haitao, Shi Yongjun, Feng Lei, Chen Xiangjun, Chen Hao, Zhang WenxiuSu Jinbo, Wang Qiong, Wang Haitao, Shi Yongjun, Feng Lei, Chen Xiangjun, Chen Hao, Zhang Wenxiu. Research on large volume airgun source signal reception in Hutubi, Xinjiang using seismic network data[J]. Earthquake Research in China, 2016, 30(3): 326332. 
Su Jinbo, Wang Qiong, Zhang Wenxiu, Wei Yunyun, Chen Hao, Wang HaitaoSu Jinbo, Wang Qiong, Zhang Wenxiu, Wei Yunyun, Chen Hao, Wang Haitao. The seasonal variation of large volume airgun signals in Hutubi, Xinjiang[J]. Earthquake Research in China, 2019, 33(2): 186194. 
Tan Junqing, Yang Runhai, Wang Bin, Xiang YaTan Junqing, Yang Runhai, Wang Bin, Xiang Ya. Weak seismic signal extraction based on the curvelet transform[J]. Earthquake Research in China, 2019, 33(2): 220234. 
Tian Xiaofeng, Yang Zhuoxin, Wang Baoshan, Yao Huajian, Wang Fuyun, Liu Baofeng, Zheng Chenglong, Gao Zhanyong, Xiong Wei, Deng XiaoguoTian Xiaofeng, Yang Zhuoxin, Wang Baoshan, Yao Huajian, Wang Fuyun, Liu Baofeng, Zheng Chenglong, Gao Zhanyong, Xiong Wei, Deng Xiaoguo. 3D seismic refraction traveltime tomography beneath the middlelower Yangtze River region[J]. Seismological Research Letters, 2018, 89(3): 9921002. DOI:10.1785/0220170245 
Wang Baoshan, Zhu Ping, Chen Yong, Niu Fenglin, Wang BinWang Baoshan, Zhu Ping, Chen Yong, Niu Fenglin, Wang Bin. Continuous subsurface velocity measurement with coda wave interferometry[J]. Journal of Geophysical Research, 2008, 113(B12): B12313. DOI:10.1029/2007JB005023 
Wang Baoshan, Ge Hongkui, Wang Bin, Wang Haitao, Zhang Yuansheng, Cai Huiteng, Chen YongWang Baoshan, Ge Hongkui, Wang Bin, Wang Haitao, Zhang Yuansheng, Cai Huiteng, Chen Yong. Practices and advances in exploring the subsurface structure and its temporal evolution with repeatable artificial sources[J]. Earthquake Research in China, 2016, 30(3): 284297. 
Wang Bin, Li Xiaobin, Liu Zifeng, et alWang Bin, Li Xiaobin, Liu Zifeng, et al. The source and observation system of Binchuan Earthquake Signal Transmitting Seismic Station and its preliminary observation results[J]. Earthquake Research in China, 2016, 30(3): 316325. 
Wang Weitao, Wang Baoshan, Jiang Shengmiao, Hu Jiupeng, Zhang YuanshengWang Weitao, Wang Baoshan, Jiang Shengmiao, Hu Jiupeng, Zhang Yuansheng. A perspective review of seismological investigation on the crust at regional scale using the active airgun source[J]. Journal of Seismological Research, 2017, 40(4): 516524 (in Chinese with English abstract). 
Wang Baoshan, Tian Xiaofeng, Zhang Yunpeng, Li Yulan, Yang Wei, Zhang Bo, Wang Weitao, Yang Jun, Li XiaobinWang Baoshan, Tian Xiaofeng, Zhang Yunpeng, Li Yulan, Yang Wei, Zhang Bo, Wang Weitao, Yang Jun, Li Xiaobin. Seismic signature of an untuned largevolume airgun array fired in a water reservoir[J]. Seismological Research Letters, 2018, 89(3): 983991. DOI:10.1785/0220180007 
Wang Weijun, Zhou Qingyun, Kou Huadong, Liu Guiping, Zhu Hongbo, Yan KunWang Weijun, Zhou Qingyun, Kou Huadong, Liu Guiping, Zhu Hongbo, Yan Kun. Improving airgun signal detection with smallaperture seismic array in Yunnan[J]. Earthquake Research in China, 2019, 33(2): 248264. 
Wei Congxin, Qin Manzhong, Zhang Yuansheng, Zou Rui, Wang Lei, Guo Xiao, Liu Xuzhou, Wang Yahong, Sun DianfengWei Congxin, Qin Manzhong, Zhang Yuansheng, Zou Rui, Wang Lei, Guo Xiao, Liu Xuzhou, Wang Yahong, Sun Dianfeng. Airgun excitation experiments at different placement depths in the Qilian Mountain of Gansu Province, China[J]. Seismological Research Letters, 2018, 89(3): 974982. DOI:10.1785/0220170253 
Wu Anxu, Ye Beng, Li Hong, Li Yun, Zheng Yiwen, Zhang Hongyan, Li YanWu Anxu, Ye Beng, Li Hong, Li Yun, Zheng Yiwen, Zhang Hongyan, Li Yan. The preliminary analysis of reliability on weak signal extraction for airgun source surveys[J]. Earthquake Research in China, 2016, 32(2): 319330 (in Chinese with English abstract). 
Yang Wei, Wang Baoshan, Liu Zhengyi, Yang Jun, Li Xiaobin, Chen YongYang Wei, Wang Baoshan, Liu Zhengyi, Yang Jun, Li Xiaobin, Chen Yong. Study on the source characteristics of downhole airgun with different excitation environments[J]. Earthquake Research in China, 2016, 30(3): 342354. 
Yang Wei, Wang Baoshan, Yuan Songyong, Ge HongkuiYang Wei, Wang Baoshan, Yuan Songyong, Ge Hongkui. Temporal variation of seismicwave velocity associated with groundwater level observed by a downhole airgun near the Xiaojiang fault zone[J]. Seismological Research Letters, 2018, 89(3): 10141022. DOI:10.1785/0220170282 
Yang Jianwen, Li Lei, Zhang Pengying, Ye Beng, He Yingwen, Chen Jia, Cha WenjianYang Jianwen, Li Lei, Zhang Pengying, Ye Beng, He Yingwen, Chen Jia, Cha Wenjian. Using airgun source signals to study regional wave velocity changes before and after the Yunlong M_{S}5.0 and Yangbi M_{S}5.1 Earthquakes[J]. Earthquake Research in China, 2019, 33(2): 320335. 
Yao Huajian, Van Der Hilst R.D., Montagner J.P.Yao Huajian, Van Der Hilst R.D., Montagner J.P. Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography[J]. Journal of Geophysical Research, 2010, 115(B12): B12307. DOI:10.1029/2009JB007142 
Zhang Yuansheng, Guo Xiao, Qin Manzhong, Liu Xuzhou, Wei Congxin, Shen Xuzhang, Yan Wenhua, Zou Rui, Yin Liang, Wang Yahong, Sun Dianfeng, Feng Hongwu, Zhang Lifeng, Guo YingxiaZhang Yuansheng, Guo Xiao, Qin Manzhong, Liu Xuzhou, Wei Congxin, Shen Xuzhang, Yan Wenhua, Zou Rui, Yin Liang, Wang Yahong, Sun Dianfeng, Feng Hongwu, Zhang Lifeng, Guo Yingxia. The construction of active source repeated monitoring in the Qilian Mountain, Gansu Province[J]. Earthquake Research in China, 2016, 30(3): 333341. 
Zheng Chenglong, Wang BaoshanZheng Chenglong, Wang Baoshan. Applications of stransform in seismic data processing[J]. Progress in Geophysics, 2015, 30(4): 15801591 (in Chinese with English abstract). 
Zhou Qingyun, Chen JunleiZhou Qingyun, Chen Junlei. Influence of different triggering conditions of airgun source on travel time changes[J]. Journal of Seismological Research, 2018, 41(2): 264272 (in Chinese with English abstract). 