Earthquake Reaearch in China  2018, Vol. 32 Issue (4): 549-559
A Study on the Seismic Velocity Changes before and after the 2016 MS6.4 Menyuan Earthquake Using the Active Source Data in the Qilian Mountain
Zou Rui, Guo Xiao, Zhang Yuansheng, Qin Manzhong, Yan Wenhua
Lanzhou Institute of Seismology, China Earthquake Administration, Lanzhou 730000, China
Abstract: The Qilian Mountain active source network data was processed using the methods of stacking, cross-correlation and interpolation, and the airgun travel time variation characteristics of P and S waves around the January 21, 2016 MS6.4 Menyua, Qinghai earthquake. The results show that about 6 months before the earthquake, the relative travel time of three stations near the epicenter showed a declined change (travel time decrease), and such a change of low value anomaly was recovered about 3 months before the earthquake. The travel time decrease then appeared again, and the earthquake occurred during the recovery process. The maximum decrease of the S-wave travel time was 18ms, and the change in travel time returned to normal after the earthquake. The variation trend of the 3 stations is consistent, including the S-wave travel time change of station ZDY38, which is nearest to the epicenter and changed obviously, and the variation range of the travel time is smaller at the stations afar. This variation pattern is related to the position of the seismic source. The shorter travel time means the velocity increase, which may be related to the regional stress accumulation.
Key words: The Menyuan, Qinghai MS6.4 earthquake     Airgun excitation signal     Travel time delay     Wave velocity variation

INTRODUCTION

With the development of observation technology and the densification of seismic stations, more studies have been conducted by seismologists on the propagation laws of natural seismic waves and the structure of the earth's interior. Earthquakes are usually considered to be caused by the accumulation and release of underground stress. Therefore, the stress state and its change of the underground medium are key factors to an earthquake (Chen Meng, 2015), and the change in the stress state of the medium will cause the velocity change of seismic waves (Birch, 1960, 1961; Scholz, 1968; Nur et al., 1969). So far, seismic waves are the only known vibrations that can penetrate the earth's interior (Chen Yong et al., 2005). Therefore, it is of great importance for earthquake prediction to understand the evolution of stress state in underground media by measuring the wave velocity changes before and after an earthquake.

At present, scholars at home and abroad use repeated earthquakes, noise and artificial sources to study the dynamic change of underground media (Schaff et al., 2004; Brenguier et al., 2008; Siliver et al., 2007; Niu et al., 2008; Chen Yong et al., 2006, 2007a; Lin Jianmin et al., 2006; Wang et al., 2008; Wang Weitao et al., 2009). However, the time and place of repeated earthquakes are uncontrollable, and the location accuracy of a natural earthquake is limited, so it is impossible to carry out long-term dynamic monitoring. The energy of the noise source is also weak and needs to be superimposed for a long time to obtain reliable measurement results, so all these factors limit the accuracy and resolution of measuring the temporal and spatial variations of underground media by natural sources (Wang Baoshan et al., 2011). While the artificial source has the advantages of known source position, controllable excitation time, flexible and convenient observation system distribution, and intensive observation, it is expected to make up for the lack of accuracy of natural sources in regional-scale research (Yang Wei et al., 2013). The use of artificial seismic source for active seismic detection has become one of the methods to detect the change of wave velocity of underground mediums with high precision, and the use of large capacity airgun source as an active source to detect and monitor the underground medium on a regional scale has the advantages of high repeatability, strong energy, remote detection distance, and also enviromental protection.

Generally, the change in the seismic wave velocity caused by stress variation of underground medium is small, so it is necessary to continuously improve the accuracy of seismic wave velocity measurements. Among several commonly used methods for measuring wave velocity, the interferometry method is the most accurate for measuring wave velocity, and the interferometric method based on cross-correlation delay detection has been widely used in seismic wave velocity measurements, and significant progress has been made. Using the active source method in the SAFOD well on the San Andreas fault in the United States of America and with the cross-correlation delay detection technique, Niu et al. (2008) successfully observed the wave velocity changes caused by atmospheric pressure and the solid tide at 1 meter depth of the well. The wave velocity increased obviously (8%) before a M3.0 earthquake occurring at 3km away from the well. In the active source detection in Kunming, Yunnan Province, Wang et al. (2008) also observed the good correlation between the variation of the medium wave velocity and atmospheric pressure variation by the coda interference method based on the cross-correlation delay detection technique.

The seismic wave travel time is closely related to the wave velocity. Accurate measurement of travel time changes can provide relevant information on the earthquake preparation process. Using the airgun excitation data of the Qilian Mountain active source from July 2015 to May 2016, selecting the airgun excitation signals received by the relevant stations in the vicinity of the MS6.4 earthquake in Menyuan, Qinghai Province in 2016, and with the cross-correlation delay detection technique, this paper measures the travel time changes and preliminarily analyzes the variation of the wave velocity before and after the Menyuan earthquake.

1 DATA SELECTION

Since the operation of the Qilian Mountain active source excitation system from July 2015, continuous excitation has been implemented every Monday evening to Tuesday early morning in order to minimize the interference of noise on the excitation signal. Each excitation interval is approximately 12 minutes, and data for 40 to 50 excitations can be obtained per week. Among them, a 24hr excitation experiment was conducted from September 30 to November 6, 2015. Up to July 1, 2016, there were a total of 6400 cumulative effective excitations, and the data is basically continuous, which provides a reliable data foundation for subsequent research.

On January 21, 2016, an MS6.4 earthquake occurred in Menyuan, Qinghai Province. The epicenter of the earthquake is located at the edge of the Qilian Mountain active source monitoring area, 178km away from the active source excitation site. The nearest active source station (station ZDY38) is 25km from the epicenter. According to the res of Bai Chaoying (1999), the source area of a magnitude 6.0 earthquake is generally tens of kilometers, therefore, we selected the stations ZDY38, ZDY37, the Himalayan second-phase mobile station No.62430 within the seismic source region, and stations ZDY28, ZDY30 and ZDY31, which are far from the epicenter (closer to the airgun source), for processing and comparative analysis (Fig. 1, Fig. 2).

 Fig. 1 Location of stations, active source and the source region of the Menyuan earthquake Red star stands for the active source excitation site, black triangle is for the active source field observation stations, and blue triangle is for the Himalayan field observation station

 Fig. 2 Travel time profiles of the selected stations by 50 times of linear superposition and 3Hz-7Hz band-pass filtering
2 PRINCIPLE AND METHOD 2.1 Principle

Moving-window cross-correlation is the most commonly used method in cross-correlation delay detection. The principle is to compare two columns of highly similar waveforms, and the reference waveform interferes with another column of waveforms. Given a time window, the maximum value of cross-correlation function is obtained by moving the time window to obtain travel time delay and wave velocity variation. The initial waveform of the recording point is defined as u(t') the waveform after disturbance is defined as $\tilde u$(t). The length of the selected time window is T. The cross-correlation function R in the time window T before and after the disturbance can be expressed as

 $R({t_s}) = \frac{{\smallint _{t - \frac{T}{2}}^{t + \frac{T}{2}}u({t^{'}})\tilde u({t^{'}} + {t_s}){\rm{d}}t}}{{\sqrt {^{'}\smallint _{t - \frac{T}{2}}^{t + \frac{T}{2}}{u^2}({t^{'}}){\rm{d}}{t^{'}}\smallint _{t - \frac{T}{2}}^{t + \frac{T}{2}}{{\tilde u}^2}({t^{'}}){\rm{d}}{t^{'}}} }}$ (1)

where, ts is the time variable; t is the central position of time window T, and when ts=t and cross-correlation function takes the maximum value, t is called the travel time delay in time window T.

When the travel time dv changes in the medium with wave velocity v, the travel time delay dt has the following relation with the elapsed time t.

 $\frac{{{\rm{d}}v}}{v} = - \frac{{{\rm{d}}t}}{t}$ (2)

Assuming that the velocity in the medium changes uniformly, and the velocity of the medium decreases (or increases), the travel time delay will increase linearly (or decrease) with the elapsed time (Snieder, 2006). The negative number of the linear slope between the linear fitting time delay and the elapsed time is the relative change of the medium's wave velocity.

2.2 The Data Processing Methods 2.2.1 Using Cross-Correlation Technique to Obtain the Moment of Excitation

The data from July 2015 to May 2016 were screened to manually remove individual waveforms with poor signal to noise ratio and recorded errors. Since some non-media changes may have an effect on the signal in the waveform (such as DC quantity, trend item, etc.), it is necessary to perform de-averaging and de-trending processing on the data at the same time, then use the cross-correlation technique to obtain the excitation time. The station ZDY22 adjacent to the Qilian Mountain active source launching site (about 80m from the airgun excitation source) is used as the reference station, and the excitation signal received by the reference station at a certain night is taken as the original template. The daily data of the reference station is cross-correlated and the excitation time is obtained (Fig. 3).

 Fig. 3 Schematic diagram of cross-correlation delay detection (a) Reference template (ZDY38). (b) Green's function (ZDY38, 2015-07-09). (c) Correlation coefficient. (d) Travel time delay before interpolation. (e) Travel time delay after interpolation. Vertical dashed line indicates the time window of the selected valid signal
2.2.2 Data Preprocessing

When the excitation time is zero, the data of each station is intercepted for 200s after the excitation time, and the sampling rate is 100Hz. In order to eliminate the influence of triggering error, water level and other factors on the airgun source, the data intercepted by each station is de-convolved separately from the excitation signal received by the reference station. Since the single-shot airgun signal will be annihilated in the noise, it is necessary to increase the signal to noise ratio by superposition to obtain an effective airgun excitation signal, so all of the data from the selected stations are superimposed as a reference template. The Green's function is obtained by superimposing the data from the stations every 2 weeks in the source area, and the Green's function with a time interval of 7 days is obtained. For the stations closer to the active source, the Green's function with a time interval of 7 days can be obtained. During the uninterrupted excitation experiment from October to November 2015, the data stacking time interval is 3 days. Finally, by bandpass filtering to 3Hz-7Hz to reduce the noise effect.

2.2.3 Cross-correlation Delay Detection

By moving the window correlation between the Green's function obtained every 2 weeks and the reference template of the corresponding station, the travel time delay can be obtained when the maximum correlation coefficient occurs. The sampling rate of the raw data is 100Hz, which means that the minimum travel time delay obtained is 0.01, but sometimes the highest peak of the correlation function is not necessarily at the sampling point, but a certain distance from the sampling point. To this end, the cosine interpolation method to restore the peak position of the correlation function is chosen, so as to obtain the time precision of the sub-sampling point (Wang Weitao et al., 2009). Fig. 3(d) and Fig. 3(e) are the travel time delay diagrams before and after interpolation. It can be seen from Fig. 3(d) and Fig. 3(e) that a higher precision travel time delay can be obtained by using cosine interpolation. According to the formula (2), the average value and the variance are obtained for the travel time delay, that is, the relative change amount and the calculation error of the relative average wave velocity are obtained.

3 MEASUREMENT RESULTS OF TRAVEL TIME BEFORE AND AFTER THE EARTHQUAKE

Station ZDY38 is located 25km away from the epicenter of Menyuan earthquake. Three time windows of 33.06s to 33.94s correspond to P-wave maximum amplitude, 56.46s to 57.60s correspond to S-wave maximum amplitude. 58.26s to 59.37s for S-wave are selected, and the travel time is calculated by cross correlation delay detection (Fig. 4). As shown in Fig. 4, the low value anomalies that lasted for nearly 2 months from early August to the end of September 2015 showed a gradual recovery. On November 4, 2015, the change in travel time showed a second decline, with a maximum range of -18ms, the abnormality of the low value continued to the MS6.4 Menyuan earthquake of January 21, 2016, and the change of travel time gradually recovered after the earthquake, and till March 25, 2016, it recovered to normal.

 Fig. 4 Travel time variation and time window selection of the selected stations in the seismic source area

For the station ZDY37, which is 41km to the epicenter of the earthquake, the travel time changes in the three time windows of 26.02s to 26.95s for the P-wave maximum amplitude, 47.10s to 48.63s for the S-wave maximum amplitude, and 49.03s to 51.87s for the S-wave are calculated (Fig. 4). It can be seen from Fig. 4 that the travel times (Figure trend of station ZDY37 and station ZDY38) are relatively consistent. First low value anomaly occurred from the beginning of August to the end of September, 2015, and then returned to normal. On November 4, 2015, the travel time once again declined. The lowest value appeared on December 28, 2016. The maximum change in travel time was -16ms. Then the travel time began to rise. During the ascent, a magnitude 6.4 earthquake occurred in Menyuan, and it returned to normal on March 14, 2016.

For the station No.62430, which is 53.8km to the epicenter of Menyuan MS6.4 earthquake, the travel time changes in the three time windows of 38.15s to 39.04s around the P-wave maximum amplitude, 63.63s to 64.64s around the maximum amplitude of the S-wave, and 65.77s to 67.15s of the S-wave are calculated. As can be seen from Fig. 5, the first low value change in travel time appeared on July 28, 2015, the low value continued until the end of September and returned to normal on September 30, 2015. On November 4, 2015, the second low value change of travel time appeared, which is different from the two stations ZDY38 and ZDY37. The travel time of station No.62430 remained to be low before the Menyuan earthquake, with the maximum change reaching -13ms. After the occurrence of the Menyuan earthquake, the lowest value of travel time appeared on January 25, 2015, and then it began to increase gradually and returned to normal on March 7, 2015.

 Fig. 5 Travel time variation and time window selection of stations far away from the epicenter

Yang Wei et al. (2010) studied the change of wave velocity before and after the Mianzhu MS5.6 earthquake. It is believed that before the earthquake, the stress in the vicinity of the fault zone gradually accumulated and the wave velocity increased, and the release of stress after the earthquake caused the wave velocity to decrease relatively. Before the MS6.4 Menyuan earthquake in 2016, the changing trend of stations ZDY28, ZDY30 and ZDY31 which are far away from the epicenter was small, and there were no obvious low value anomalies (Fig. 5). About 6 months before the earthquake, the relative travel time of stations such as ZDY38, ZDY37, and 62430 near the epicenter showed a downward change (reduction in travel time), and the low value anomaly was recovered about 3 months before the earthquake, then the travel time decreased again. This anomaly may be related to the continuous accumulation of the stress of the underground medium in the seismogenic zone, which leads to an increase of velocity of seismic waves propagating through the region.

4 CONCLUSION AND DISCUSSION 4.1 Error Analysis

The lower error limit of the travel time delay στ can be calculated by the Cramer-Rao Lower Bound rule (Liu Zifeng et al., 2015), i.e.

 ${\sigma _\tau } \ge \sqrt {\frac{3}{{2{f_0}^3{\pi ^2}T({B^3} + 12B)}}\left[ {\frac{1}{{{\rho ^2}}}{{\left({1 + \frac{1}{{SNR}}} \right)}^2} - 1} \right]}$ (3)

Where, f0 is the main frequency signal, T the length of time window; B the ratio of signal to frequency, ρ the correlation coefficient of waveform, and SNR the ratio of signal to noise. Formula (3) shows that the similarity of the waveform, the signal to noise ratio and the ratio of the frequency width of the signal to the main frequency are the main factors that affect the measurement accuracy of the travel time delay (Zhang Jinchuan et al., 2014).

The primary frequency of the Qilian Mountain active source is basically 5Hz, and the energy of the signal is mainly concentrated at 4.5Hz to 6.5Hz. Therefore, f0≈5Hz is selected, B≈0.4. When the cross-correlation delay detection is performed, the length of the window T=0.4s, ρ≈1, SNR≈15 is selected, and the formula (3) can be simplified as

 ${\sigma _\tau } \ge \frac{1}{{2\pi {f_0} \cdot SNR}}\sqrt {\frac{1}{{{f_0}TB}}}$ (4)

From formula (4), we get the lower limit of the theoretical error of the travel time delay, which is about 2.37×10-4s. When the cross-correlation is performed, the mean square error of the travel time delay is about 6.15×10-4s. It can be seen that the measured error is very close to the theoretical error, but it is one order of magnitude smaller than the observed travel time variation value of stations ZDY38, ZDY37, and 62430.

In addition to the calculation error, the accuracy of the measurement is also affected by the observing system and environmental factors (Yang Wei et al., 2010). The Qilian Mountain active source stations adopt the GPS continuous timing method, and its timing accuracy is 10-7μs, which is three orders of magnitude higher than the calculated accuracy of the direct S-wave travel time variation value. The influence of solid tide and air pressure is a complicated problem. The superposition of data has a certain weakening effect on these influences. Therefore, the waveform data after superimposition processing may be affected by the influence of solid tide and air pressure.

4.2 Conclusion and Discussion

In this paper, we use the airgun excitation data from the Qilian Mountain active source from July 2015 to May 2016 to select the airgun excitation signals received by the relevant stations near to and far away from the epicenter of the Menyuan, Qinghai MS6.4 earthquake. The time-delay detection technique is applied to measure the travel time before and after the MS6.4 Menyuan earthquake. The results show that the trend of travel time change of the stations ZDY28, ZDY30, ZDY31 and other stations far away from the epicenter is small, and there is no obvious low-value anomaly process. About 6 months before the earthquake, the relative travel time of the three stations near the epicenter showed a declined change (the travel time decreased), and about 3 months before the earthquake, the low value anomaly recovered, and then the travel time decreased again. The MS6.4 earthquake occurred in the process of travel time recovery. The maximum decrease of S-wave travel time is -18ms, and the change of travel time gradually returned to normal after the earthquake, and the change trend of the 3 stations is more consistent.

Before the earthquake occurred, the stress in the vicinity of the fault zone gradually accumulated and the wave velocity increased, and the stress release after the earthquake caused the wave velocity to decrease relatively (Yang Wei et al., 2010). Shorter travel time means increased speed, which may have a certain relationship with regional stress accumulation. The travel time change of the distant stations to the epicenter such as ZDY28, ZDY30, and ZDY31 is small. The travel time change of the three stations ZDY38, ZDY37 and 62430 near the epicenter presents a low value before the earthquake, which may be related to the stress accumulation in the preparation process of the Menyuan MS6.4 earthquake of 2016, which may increase the velocity of seismic waves propagating through the region. The relative change of the average S-wave velocity before the earthquake is 0.38% in station ZDY38, 0.27% in the station ZDY37, and 0.15% in the station No.62430. It is clear that the station ZDY38, which is nearest to the epicenter of Menyuan earthquake and the Changma-Rubo fault, sees the most obvious change. The research on the processing and application of the observation data of the airgun active source in the Qilian Mountain has just started. The results obtained in this paper are preliminary, but the results will be helpful to the rapid production of observation data and the analysis of the information related to the change in seismic wave velocities, so as to provide service for local earthquake prevention and disaster reduction as soon as possible.

ACKNOWLEDGEMENTS

The authors would like to express their heartfelt thanks to research professor Wang Baoshan from the Key Laboratory of Seismic Observation and Geophysical Imaging of Institute of Geophysics, China Earthquake Administration, who provided the cross-correlation delay calculation program.

REFERENCES
 Bai Chaoying. Study on the relation between the moderate earthquake concentrating regions and the three parameters of future strong earthquakes[J]. Northwestern Seismological Journal, 1999, 21(1): 48–54 (in Chinese with English abstract). Birth 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. Birth 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., Nadeau R.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 Yong, Zhu Rixiang. Proposed project of "Underground Bright Lump"[J]. Advances in Earth Science, 2005, 20(5): 485–489 (in Chinese with English abstract). Chen Yong, Zhou Huawei, Ge Hongkui. Seismic array in North China[J]. Journal of Geodesy and Geodynamics, 2006, 25(4): 1–5 (in Chinese with English abstract). Chen Yong, Wang Baoshan, Ge Hongkui, Xu Ping, Zhang Wei. Proposed of transmitted seismic stations[J]. Advances in Earth Science, 2007a, 22(5): 441–446 (in Chinese with English abstract). Lin Jianmin, Wang Baoshan, Ge Hongkui, Chen Qifu, Chen Yong. Doublet and its potential application in active exploration[J]. Earthquake Research in China, 2006, 22(1): 1–9 (in Chinese with English abstract). Liu Zifeng, Su Youjin, Wang Baoshan, Wang Bin, Yang Jun, Li Xiaobin. Study on analysis method of travel time variations of seismic wave of active source in Binchuan[J]. Journal of Seismological Research, 2015, 38(4): 591–597 (in Chinese with English abstract). 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. Nur A., Simmons G. The effect of saturation on velocity in low porosity rocks[J]. Earth and Planetary Science Letters, 1969, 7(2): 183–193. DOI:10.1016/0012-821X(69)90035-1. Schaff D.P., Beroza G. Coseismic and postseismic velocity changes measured by repeating earthquakes[J]. Journal of Geophysical Research, 2004, 109(B10): B10302. DOI:10.1029/2004JB003011. Scholz C.H. Micro fracturing and the inelastic deformation of rock in compression[J]. Journal of Geophysical Research, 1968, 73(4): 1417–1432. DOI:10.1029/JB073i004p01417. Silver P.G., Daley T.M., Niu Fenglin, Majer E.L. Active source monitoring of cross-well seismic travel time for stress-induced changes[J]. Bulletin of the Seismological Society of America, 2007, 97(1B): 281–293. DOI:10.1785/0120060120. Snieder R. The theory of coda wave interferometry[J]. Pure and Applied Geophysics, 2006, 163(2/3): 455–473. Wang Baoshan, Zhu Ping, Chen Yong, Niu Fenglin, Wang Bin. Continuous subsurface velocity measurement with coda wave interferometry[J]. Journal of Geophysical Research, 2008, 113: B12313. DOI:10.1029/2007JB005023. Wang Baoshan, Wang Weitao, Ge Hongkui, Xu Ping, Wang Bin. Monitoring subsurface changes with active sources[J]. Advances in Earth Science, 2011, 26(3): 249–256 (in Chinese with English abstract). Wang Weitao, Wang Baoshan, Ge Hongkui, Chen Yong, Yuan Songyong, Yang Wei, Li Yijin. Using active source to monitor velocity variation in shallow sediment caused by the Wenchuan Earthquake[J]. Earthquake Research in China, 2009, 25(3): 223–233 (in Chinese with English abstract). Yang Wei, Ge Hongkui, Wang Baoshan, Yuan Songyong, Song Lili, Jia Yuhua, Li Yijin. Velocity changes observed by the precisely controlled active source for the Mianzhu MS5.6 earthquake[J]. Chinese Journal of Geophysics, 2010, 53(5): 1149–1157 (in Chinese with English abstract). Yang Wei, Wang Baoshan, Ge Hongkui, Wang Weitao, Chen Yu. The active monitoring system with large volume airgun source and experiment[J]. Earthquake Research in China, 2013, 29(4): 399–410 (in Chinese with English abstract). Zhang Jinchuan, Wang Qincai, Xue Bing, Ding Lisha. Monitoring media velocity variations with coda wave interferometry[J]. Earthquake, 2014, 34(3): 62–73 (in Chinese with English abstract).