Earthquake Research in China  2019, Vol. 33 Issue (4): 557-572     DOI: 10.19743/j.cnki.0891-4176.201904006
Study on Azimuth Variation of Airgun Signal Propagation
HE Suge, ZHOU Qingyun
Yunnan Earthquake Agency, Kunming 650041, China
Abstract: Based on the data recorded by the observation network during the intensive excitation period from November to December 2015 at Binchuan Earthquake Signal Transmitting Seismic Station (BESTSS) in Yunnan Province,the noise in waveform recording is removed by S-transform template filtering method,and the azimuth of airgun signal propagation is calculated and analyzed from the horizontal waveform recordings. The results show that:① the azimuth angle of airgun signal after propagation is sensitive to stress change,and can clearly reflect the diurnal variation of tidal stress,which can be used to monitor the change of stress state in crustal medium; ② the azimuth angle of airgun signal in some stations has changed abruptly after propagation on December 4,which may be related to the change of airgun source; ③ five-shot superposition or five-point smoothing of azimuth angle of single shot are carried out for airgun signals in stations far away from epicenter,and results show that azimuth angle from superposition or smoothing is more stable and has clear diurnal variation characteristics after propagation.
Key words: Airgun source     Azimuth change     S-transform     Template matched filtering     Crustal medium change

INTRODUCTION

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, signal-to-noise ratio (SNR) of seismic signals is improved by means of encrypted observation array and long-time 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 signal-to-noise, 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 high-quality 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 signal-to-noise 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 P-wave 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 P-wave velocity, S-wave 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 MS5.0 and Yangbi MS5.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 DATA

The 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, 0-12 oclock every hour per shot, 12-24 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 signal-to-noise ratio of each waveform. Therefore, we use the data of stations with less noise and clear signal of each shot.

 Fig. 1 Station distribution map Red Pentagon is the position of BESTSS. Solid triangles are stations used in this paper. Red triangles are stations with the azimuth angle increasing on December 4, greens are reduced station and blacks are the station with no obvious changes

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).

 Fig. 2 Vertical component correlation result (solid circle) and amplitude (open triangle) of the airgun signal recorded by CKT1. The epicenter distance of the station is about 40m
2 METHOD AND DATA PROCESSING 2.1 S-transform Template Filtering

Stockwell R. G. et al., (1996) proposed the S-transform method, which is a time-frequency analysis method. This method combines the advantages of short-time Fourier transform and wavelet transform. Based on S-transform, scholars have proposed different filtering methods. Zheng Chenglong et al. (2015) proposed S-transform 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 signal-to-noise ratio. The time-frequency spectrum S(τ, f) of the superimposed signal h(t) is obtained by S-transform, and the time-frequency 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 S-transform, multiplication with filter template and S-inverse transformation into time domain.

2.2 Calculation of Back Azimuth of Airgun Signal

The P-wave 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 P-wave 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 P-wave window length, and ui(t) and uj(t) are the horizontal components of the seismic record. In the absence of noise, the covariance matrix c has a non-zero 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 non-zero eigenvalues, and the ratio of the two non-zero eigenvalues reflects the linearity of the particle motion. The S-transformation template filter can significantly suppress the noise in the airgun signal, and the signal-to-noise ratio can be used to estimate the error in the azimuth calculation (Drew J. et al., 2008).

2.3 Data Processing Flow

First, preprocessing. Waveforms taken from continuous waveform recordings typically contain non-zero 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; de-averaging, removing the non-zero 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 1-10Hz, which removes some common noises and ensures that the airgun signal is not filtered.

Then, filtering using S-transform template filter. The airgun signal recorded by stations located 10km away from the epicenter is relatively low in signal-to-noise ratio and cannot be used directly. S-transform template filtering can suppress most noise under the premise of retaining effective airgun signal, and effectively improve the signal-to-noise 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 S-transform to the time domain. The filter template threshold is 0.15. It can be seen from the filtering results that the S-transform template filtering can greatly suppress the noise while effectively retaining the target signal (Fig. 3).

 Fig. 3 An example of S-transform template filtering (a) Single shot signal of Station 53257. (b) Superposition signal of Station 53257. (c) Filter template. (d) Single shot signal filtered by S-transform template

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 P-wave arrives at the station (P-wave 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 near-source stations and the arrival time of the S-wave at the far-reaching stations. Therefore, we use the termination time as follows:

 ${t_{{\rm{end }}}} = {t_{{\rm{porrive }}}} + {\rm{dist}}/14 + 3.25$ (3)

Where tp_arrive is the time when the P-wave arrives at the station, and dist is the epicenter distance of the station.

3 DATA PROCESSING RESULTS

When 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 signal-to-noise 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:

 Fig. 4 Changes of back azimuth angle of Station 53258 and other stations

 Fig. 5 Changes of back azimuth angle of Station 53254 and other stations

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.

 Fig. 6 Variation of back azimuth of some stations from November 18 to 28

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 Analysis

In the absence of noise, the azimuth calculated by equation (2) is the direction of P-wave vibration; when there is noise, the calculated azimuth may not coincide with the direction of P-wave 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 100.1 to 104. We found that the higher the signal-to-noise 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).

 Fig. 7 The influence test of signal-to-noise ratio on azimuth calculation of EY02 station The added noise is 1646 segments of noise and 40 gain levels

After S-transform template filtering, the amplitude of the noise portion is set to 0, so the traditional method cannot be used to calculate the signal-to-noise 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 P-wave arrival, filter the noise, and use the filtered data to evaluate the noise level; after P-wave 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 Sj is the value of the jth sampling point of the airgun signal, and Nni is the value of the ith sampling point of the nth segment of noise data.

We calculated the signal-to-noise 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 signal-to-noise 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 signal-to-noise 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.

 Fig. 8 (a) Signal-to-noise ratio (SNR) of every shot at each station in Fig. 3(b). (b)Signal-to-noise ratio (SNR) of every shot at each station in Fig. 4
4.2 Analysis of the Cause of Azimuth Variation

There are pores in the rock. When the intensity or direction of the stress in the rock changes, the open-close degree and direction of the rock pores will change (Crampin S., 1987), in which the change of pores open-close 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 1-3 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 P-wave 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 v1 and v2 are the P-wave 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 v1=3.0km/s, v2=4.5km/s, α1=1.2×10－6Pa－1, α2=4×10－7Pa－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 P-wave 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 signal-to-noise 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 signal-to-noise ratio is low; human activities are weakened at night, reduced noise intensity, and increased signal-to-noise 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 S-transform 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 signal-to-noise 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 P-wave 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 4

The 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.

 Fig. 9 Azimuth changes of stations and airgun signal correlation coefficient from December 4 to 5, 2015
4.4 Data Processing Method for Long Distance Stations

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 single-shot 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 S-transform template filtering method is used to filter the single-shot 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.

 Fig. 10 The calculation results of back azimuth angle obtained by different processing methods of Stations 53251 and 53252 (a) The results of single-shot calculation at 53251 station. (b) The results of 5-point smoothing at 53251 station. (c) The results of 5-shot superposition at 53251 station. (d) The results of single-shot calculation at 53251 station. (e) The results of 5-point smoothing at 53252 station. (f) The results of 5-shot superposition at 53252 station
4.5 Threshold Selection of Filtering Template

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.

 Fig. 11 The correlation coefficients between the calculated results of different filtering templates Δ is the correlation coefficient between the filtering result of threshold value of 0.05 and the filtering result of 0.15, and ○ is the correlation coefficient between the filtering result of threshold value of 0.15 and the filtering result of 0.30
4.6 Reasons for the Difference between Calculated Azimuth and Real Azimuth

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 P-wave phase should be the best selected. Since the source time function of the Binchuan active source signal lasts longer than 3s, the P-wave and S-wave 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 low-speed rock (such as clastic rock Qp, mudstone J3) and hard high-speed rock (such as limestone D1, 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 long-distance 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 S-transform Template Filtering on Azimuth Calculation

Filtering, 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 S-transform template filtering and non-S-transform 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 S-transform template filtering method can improve the waveform quality and stability, and has little influence on the azimuth variation analysis.

 Fig. 12 Comparison between the results of S-template filtering with those of non-S-template filtering (a) Azimuth results of 53273 stations. (b) Correlation coefficients between single-shot waveforms and superimposed waveforms of 53273 stations. (c) Azimuth results of EY01 stations. (d) Correlation coefficients between single-shot waveforms and superimposed waveforms of EY01 stations
5 CONCLUSION

(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 moderate-strong 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 S-transform template filtered single-shot signal is poor. The 5-shot superposition or 5-point smoothing method can effectively improve the stability of the calculated results.

REFERENCES
 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): 298-315. Birch F.Birch F. The velocity of compressional waves in rocks to 10 kilobars:1[J]. Journal of Geophysical Research, 1960, 65(4): 1083-1102. 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): 2199-2224. 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): 1478-1481. 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 large-volume airgun source excitation conditions in a reservoir[J]. Earthquake Research in China, 2016, 30(3): 355-363. 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): 281-283. 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): 1739-1751. DOI:10.1007/s11430-016-9096-6 Cheng Xin, Niu Fenglin, Wang BaoshanCheng Xin, Niu Fenglin, Wang Baoshan. Coseismic velocity change in the rupture zone of the 2008 MW7.9 Wenchuan earthquake observed from ambient seismic noise[J]. Bulletin of the Seismological Society of America, 2010, 100(5B): 2539-2550. 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): 6147-6164. DOI:10.1029/JZ070i024p06147 Crampin S.Crampin S. Geological and industrial implications of extensive-dilatancy anisotropy[J]. Nature, 1987, 328(6130): 491-496. 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: 1446-1450. 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): 125-135. DOI:10.1111/j.1365-246X.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): WC143-WC155. DOI:10.1190/geo2010-0409.1 Hillers G., Ben-Zion Y., Campillo M., Zigone D.Hillers G., Ben-Zion 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): 920-932. 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): 235-247. 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): 305-319. 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): 511-520. 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): 36-45. 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): 204-208. DOI:10.1038/nature07111 Niu Fenglin, Li JuanNiu Fenglin, Li Juan. Component azimuths of the CEArray stations estimated from P-wave particle motion[J]. Earthquake Science, 2011, 24(1): 3-13. 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 MW7.9 2008 Wenchuan earthquake[J]. Journal of Geophysical Research, 2014, 119(4): 3155-3168. 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): 5719-5731. 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 middle-lower Yangtze river region in eastern China from surface-wave tomography of a large volume airgun-shot experiment[J]. Seismological Research Letters, 2018, 89(3): 1003-1013. 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): 174-185. 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): 998-1001. 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): 326-332. 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): 186-194. 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): 220-234. 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 travel-time tomography beneath the middle-lower Yangtze River region[J]. Seismological Research Letters, 2018, 89(3): 992-1002. 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): 284-297. 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): 316-325. 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): 516-524 (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 large-volume airgun array fired in a water reservoir[J]. Seismological Research Letters, 2018, 89(3): 983-991. 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 small-aperture seismic array in Yunnan[J]. Earthquake Research in China, 2019, 33(2): 248-264. 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): 974-982. 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 air-gun source surveys[J]. Earthquake Research in China, 2016, 32(2): 319-330 (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): 342-354. Yang Wei, Wang Baoshan, Yuan Songyong, Ge HongkuiYang Wei, Wang Baoshan, Yuan Songyong, Ge Hongkui. Temporal variation of seismic-wave velocity associated with groundwater level observed by a downhole airgun near the Xiaojiang fault zone[J]. Seismological Research Letters, 2018, 89(3): 1014-1022. 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 MS5.0 and Yangbi MS5.1 Earthquakes[J]. Earthquake Research in China, 2019, 33(2): 320-335. 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): 333-341. Zheng Chenglong, Wang BaoshanZheng Chenglong, Wang Baoshan. Applications of stransform in seismic data processing[J]. Progress in Geophysics, 2015, 30(4): 1580-1591 (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): 264-272 (in Chinese with English abstract).