Earthquake Research in China  2019, Vol. 33 Issue (4): 584-595     DOI: 10.19743/j.cnki.0891-4176.201904009
Comparative Study on Two Methods for Measuring Wave Velocity Change of Crustal Medium Based on Large Volume Airgun Excitation Data
YE Beng1,2, CAO Wenzhong2, HUO Yuanhang2, CHEN Jia1, LI Xiaobin1
1. Office of the Western Yunnan Earthquake Prediction Study Area, CEA, Dali 671000, Yunnan, China;
2. University of Science and Technology of China, Hefei 230000, China
Abstract: This paper proposes the application of dynamic programming method to calculate the relative change of wave velocities and compares its similarities and differences with the cross-correlation delay estimation method based on interference. The results show that:① the trend of wave velocities obtained by cross-correlation method and dynamic programming method are consistent. Besides,it is considered that the calculated result using cross-correlation delay method is reliable. ② Compared with the cross-correlation delay method,the calculated result of the dynamic programming method has a magnifying effect and is more sensitive to small disturbances. ③ Under ideal conditions,the wave velocity change trend calculated by P-wave and S-wave phase should be consistent. In addition,the cross-correlation delay method is used to calculate the wave velocity change. Under appropriate conditions,the process of recovering from the suspected wave velocity before the ML1.1 earthquake near the airgun source can be observed.
Key words: Artificial source     Airgun source     Cross-correlation     Wave velocity change     Dynamic programming

INTRODUCTION

Earthquakes are destructive natural disasters for human beings. Every time a strong earthquake occurs, it will cause huge losses. It is our urgent hope to understand the origin and occurrence of earthquakes and thus reduce the loss of earthquake disasters. Due to the complex physical processes of earthquakes, the inaccessibility of the earths interior and other factors hinder the understanding of the earthquake-triggered environment, and finding a physical quantity with a clear meaning related to earthquake is the primary condition for predicting earthquakes (Wang Weitao et al., 2009). Wave velocity is jointly determined by the intrinsic properties of a rock and the external environment in which it is located. The mineral composition and content of a rock, physical properties, distribution state, porosity, pore distribution state, fluid properties in the pores, temperature, pressure, and pore pressure all have great influence on the change of wave velocity (Yang Xiaosong et al., 2003), which is a suitable physical quantity for describing the change of the properties of the medium in the seismogenic environment. The change in wave velocity reflects the change in the physical properties of the subsurface medium. Therefore, by measuring the change of wave velocity, the spatial and temporal changes of underground media can be inferred, and their relationships with the occurrences of earthquakes are studied.

The predecessors have conducted extensive researches on the wave velocity changes before and after earthquake seismic event. For example, after Japanese scholars used the interferometry method to measure the effects of seasonal factors, the surface shear wave velocity varies 5% before and after the 2011 "3·11" earthquake in Japan (Nakata N. et al., 2011). Numerous studies have been done in this area by Chinese scholars. Xu Z. J. et al., (2009) found that after three earthquakes in Sumatra in 2004, 2005 and 2007, the Rayleigh wave velocity in the Sumatra region changed significantly. Additionally, Cheng Xin (2010) used noise-related method to measure the wave velocity of the Rayleigh wave on the northwest side of the Longmenshan fault after the Wenchuan earthquake and found it decreased by 0.4%. Liu Mian et al. (2014) further confirmed that the Longmenshan fault had a corresponding change in stress before the Lushan earthquake. Pei Shunping et al. (2019) also indicated that the body wave velocity in the area decreased to some extent before the Wenchuan earthquake by picking up the body waves of small earthquakes. The difficulty in studying the wave velocity change before and after large earthquakes is the accuracy of the wave velocity measurement. Previously, the reliable high-precision data of wave velocity was obtained from rock physics experiments, but the measurement of high-precision wave velocity variation at the regional scale has always been a difficult point in geophysics (Lin Jianmin et al., 2006). Since the events generated by active source have the advantages of position controllability, high repeatability, and large detection scale (Wang Bin et al., 2016; Yang Wei et al., 2013), while the natural seismic positioning accuracy is limited, and the noise source energy is weak, the former one is utilized oftentimes. The wave velocity measurement method is able to enhance the measurement accuracy, thus has attracted the attention of many scholars. In the early periods, the active sources that geophysicists used were airguns (Yang Wei et al., 2013), ultrasound (Wang Zijie et al., 1997), piezoelectric ceramics (Niu Fenglin et al., 2008), repeated earthquakes (Zhou Longquan et al., 2007), and repeated blasting (Li Le et al., 2007). However, with the development of the active source technology and high-precision observation systems, airgun has received much attention as a important source so that large-capacity airgun excitation sources have been widely used in China since the 21st century.

The variation of wave velocity is extracted from the experimental data collected by airgun source. It is an important scientific issue to study the change of crustal medium before and after an earthquake (Chen Yong et al., 2007). Studies suggested that there are abnormal fluctuations in wave velocity before and after an earthquake, and it may related to factors such as solid tide, atmospheric pressure, and rainfall (Silver P.G. et al., 2007; Wang Baoshan et al., 2008).

The regional velocity fluctuation of the subsurface medium is extracted from the seismic waveform obtained by airgun active source using the cross-correlation delay method (Wang Weitao et al., 2009; Luo Guichun, 2006; Wang Bin, 2009). Compared with previous studies, this method makes use of the advantages of the airgun source system, i.e. the high repeatability of the source and receiving stations, making themeasurement more straight forward. In practical applications, it is necessary to set a calculation time window, and it is assumed that the amplitude of the phase in the time window is constant for different events. Since the waveform component of the active source signal generated by airgun is complex, when the displacement amount is large and the amplitude changes drastically, the method may fail (Liu Lanfeng et al., 2015; Hale D., 2013; Hall S. A., 2006). This paper introduces the dynamic programming methods, and find which method is more advantageous by compares the calculation results of the two methods.

1 METHODOLOGY 1.1 Cross-Correlation Delay Method

Due to the great uncertainty of the occurrence of natural earthquakes, there are numerous problems in directly measuring the absolute wave velocity changes. However, in recent years, it has become possible to use artificial sources to measure wave velocity changes. Among them, Niu Fenglin et al. (2008) is the most famous one which implemented the measurement work in the St.Anders fault deep exploration program. Both the airgun source excitation environment and the obtained phase are complex. Generally, when calculating the wave velocity, a cross-correlation-based measurement technique is used, that is, measuring the delay time of a certain phase multiple times in a fixed path, and measuring accuracy can be higher than the sub-sampling rate and achieve the desired accuracy (Wang Baoshan et al., 2011). The principle is as follows:

For the active source excitation signal repeated twice in the same place, the waveform information correlation function is expressed as

 ${R_{xy}}(\tau) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} x (t){y^*}(t - \tau){\rm{d}}t$ (1)
 ${R_{yx}}(\tau) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} y (t){x^*}(t - \tau){\rm{d}}t$ (2)

Among them, the seismic waveform record is often a power-limited signal, and the two pieces of information are represented as the sum of the two real numbers, x(t) and y(t), and there is a delay τS between them.

 ${R_{xy}}\left({{\tau _{\rm{S}}}} \right) = \frac{{\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} x (t)y(t - \tau){\rm{d}}t}}{{\sqrt {\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} {{x^2}} (t){\rm{d}}t\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} {{y^2}} (t){\rm{d}}t} }}$ (3)

For the airgun active source signal received twice by the same station, the real waveform x(t) and the waveform y(t) are highly similar, but there is time delay, that is, when the correlation coefficient takes the maximum value, the corresponding τS is the delay of the two signals.

 ${R_{xy}}\left({{\tau _{\rm{S}}}} \right) = \frac{{\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} x (t)y\left({t - {\tau _{\rm{S}}}} \right){\rm{d}}t}}{{\sqrt {\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} {{x^2}} (t){\rm{d}}t\int_{t - \frac{T}{2}}^{t + \frac{T}{2}} {{y^2}} (t){\rm{d}}t} }}$ (4)

In the analysis, the experimental signal is selected as the standard value, and the delay between the remaining excitation signal and the standard signal indicates the relative travel time difference between the two excitation events. For the same phase, the path from the active source to the station is fixed, and the travel time difference represents the speed variation of the phase as it passes through the medium.

1.2 Dynamic Programming

The dynamic specification method was first applied to the field of speech recognition (Sakoe H. et al., 1978) to identify waveform changes caused by the pronunciation of the same words. Anderson K. R. et al., (1983) first implemented it in the field of high-precision seismic imaging research. Anderson discussed the global path constraints on the basis of predecessors, and proposed that higher slope accuracy can be obtained when the slope constraint is reduced. Keysers D. et al., (2013) introduced it into 2D seismic imaging and found that the method is computationally intensive. Pishchulin L., (2010) proposed an approximate calculation that reduces the amount of computation. Based on Mottls approximation method, Hall S. A., (2006) proposed to divide the optimization problem into two steps: accumulation and backtracking. He also gave the mathematical solution process and deducted it into multidimensional image calculation. Liu Lanfeng et al. (2015) used Daves method to calculate the multidimensional space and found that dynamic image deformation is less sensitive to noise than local correlation imaging. In this theory, the dynamic programming method finds a path through a number of grid points, and the cumulative minimum value of the grid points passing through the path is the optimal value. In the calculation steps, only the matching distance of the frame is calculated. It satisfies the time correspondence of the time warping function description template and the reference template as a constraint to solve the regular function corresponding to the minimum cumulative distance when the two templates match. The process is listed below:

Suppose there are two series f[i] and g[i].

 $f=f_{1}, f_{2}, \cdots \cdots, f_{i}$ (5)
 $g=g_{1}, g_{2}, \cdots \cdots, g_{j}$ (6)

To find the similarity of two series, we can build a matrix, i.e.

 $c(k)=(f(k), g(k))$ (7)

The above formula represents the pairing of any two points in the two series.

 $\mathrm{d}(c(k))=\mathrm{d}(f(k), g(k))=\mathrm{d}(i, j)=\left\|f_{i}-g_{j}\right\|$ (8)

Wherein, the matrix element d(c(k)) represents the distance between any two points in the two series.

 $D(f, g) = \frac{1}{N}\min F\left[ {\sum\limits_{k = 1}^k {\rm{d}} (c(k)) \cdot w(k)} \right]$ (9)

where, D(f, g) is all accumulation. When it reaches the minimum, it means that the two series have the maximum similarity. w(k) is the weight of each step. Here, w(k)=1 and it means all the steps.

According to the theory, the constraints are:

(1) Relationship between the previous step and the latter step

 \begin{aligned} i(k-1) \leqslant i(k) & \text { and } j(k-1) \leqslant j(k) \\ \end{aligned} (10)
 \begin{aligned} i(k)-i(k-1) & \leqslant 1 \text { and } j(k)-j(k-1) \leqslant 1 \end{aligned} (11)

that is

 $c(k-1)=\min \left\{\begin{array}{c} {(i(k), j(k)-1)} \\ {(i(k)-1, j(k)-1)} \\ {(i(k)-1, j(k))} \end{array}\right\}$ (12)

(2) Boundary conditions

 $i(1)=1, j(1)=1, i(k)=I, j(k)=J$ (13)

(3) Sliding window size before and after

 $|i(k)-j(k)| \leqslant r$ (14)

Where i and j represent the lengths of f[i] and g[i], respectively. r is the sliding window set according to the actual situation.

Through the above conditions, the change rate of the time shift amount is constrained to find a global optimal time shift amount. That is, when the minimum value D(f, g) is recorded, the relative time difference of the corresponding two series windows is recorded, and the time delay can be obtained.

In this paper, since seismic imaging is not required, the amount of data calculation is relatively less and the approximate calculation method is not used.

As the theory shown above, the dynamic programming method can tolerate a greater degree of compression or tensile deformation of the data when the correlation coefficient is obtained.

2 TEMPLATE INSPECTION 2.1 Building a Theoretical Model

There are many seismic wave simulation methods in complex media (Zhang Wei, 2006). This paper uses a purely numerical finite difference method based on lattice points. The data is generated by two-dimensional forward modeling to perform a test on the two wave velocity variation calculations above. The relevant parameters of the media model are set as follows.

According to the parameters listed in Table 1, a simple medium model of the Binchuan area is constructed. The Rake wavelet is set as the source on the surface of the model, and 10 receivers are arranged with an interval of 1km at the same level as the source. A medium-disturbing body is set 5km from the source. The body is a rectangle which is 3km long and 1km wide that moves between 2km and 23km in depth. The disturbing body moves down 1km each time during the experiment. The media model and the received data are shown in Fig. 2.

Table 1 Theoretical model related parameters

 Fig. 1 (a) Corresponding point-point diagram of the cross-correlation calculation of similarity time.(b) Corresponding point-point diagram of the dynamic programming calculation of similarity time

 Fig. 2 Waveforms recorded by the media model and receiver (a) Media model related parameters and receiving station distribution. (b) All event waveforms received by Station 3;All event waveforms received by Station 9;Waveforms received by the disturbing bodies at all stations with a depth of 3km

After 24 simulation experiments, the data from the stations with the epicenter distances of 6km and 12km are selected and analyzed. The results show that during the sinking process of the disturbing body, multiple waves are obviously deformed, while body waves are not. For body waves, the influence of the disturbing body on S-wave is greater than that on P-wave (Fig. 2).

2.2 Analysis of Experimental Results

The cross-correlation delay method and the dynamic programming method are used to calculate the relative change of the wave velocity. However, in the two calculations, the time window for calculating the relative change is unchanged. According to the ray theory, the ray paths of different phases are different. Therefore, time-lapse windows of the most easily identifiable S and P are used to calculate the wave velocity variation, that is, the red and green marker bands in the first row of Fig. 3, which shows the epicentral distances and the wave velocity variations of the stations 6km and 12km from the epicenter. The second row is the calculation result using the cross-correlation delay method, while the third row is the calculation result using the dynamic programming method. The abscissa is time and the ordinate is the event number.

 Fig. 3 The relative travel time of the P phase (red) and S phase (green) and their respective template phases in each event calculated by the theoretical waveforms (a)Template waveform for all events of station. (b) The template waveform for all events of station. (c) and (d) are the calculation results for cross-correlation delay method. (e) and (f) are the calculation results for dynamic programming method. (c) and (e)reveal the relative travel time difference between the corresponding phases of the stations with the epicentral distance of 6km and the template phase. (d) and (f) are the relative travel time difference between the corresponding phases of the stations at the epicentral distance of 12km and the template phase

As shown in Fig. 3, the cross-correlation delay method and the dynamic programming method have the following similarities and differences:

(1) The calculation result using the dynamic programming method is one order of magnitude larger than the result using the cross-correlation delay method, but the basic trend of the wave velocity change is consistent. Since the disturbing body is a rectangle which is 3km long and 1km wide, the wave velocity in the disturbing body is 9km. Therefore, the maximum time difference of two time events caused by the disturbing body is 0.11s, and the maximum time difference of the wave is 0.061s. The calculation results of the dynamic programming method are too large, while the results of the cross-correlation method based on interference are too small. In addition, the calculation results of the P-wave are close to the theoretical values. The reason for this difference may be related to their related principles, and may also be related to the paths of body waves in different events.

(2) Compared with the cross-correlation delay method, the dynamic programming method has a more obvious amplification effect and is more sensitive to small disturbances. Correspondingly, the cross-correlation method has strong stability.

(3) The wave calculation results obtained by the two methods show that the change trend of P-wave calculation is almost the same and the calculation result difference of the S-wave is relatively large but with a consistent basic trend. Besides, the calculation result of the wave is more sensitive to the disturbance. The calculation result of S-wave is more sensitive to disturbance and there is an amplification effect.

(4) The difference between the two methods increases with the epicentral distance.

In summary, from the experimental results, the cross-correlation method based on interference is closer to the theoretical value. This may be related to the amplification effect of the small changes in the waveform caused by the dynamic programming method when calculating the correlation coefficient due to its tolerance to the waveform shape change.

3 CALCULATION RESULTS AND DISCUSSION OF EXPERIMENTAL DATA IN BINCHUAN AREA

During December 12-16, 2015, the Binchuan Seismic Signal Transmitting System has continuously conducted a 4-day experiment at a frequency of one excitation per 1hour. In this period, three small earthquakes occurred in the inner center of the seismic signal launching platform. Among them, the earthquake with ML1.1 occurred at 20:46 on December 12, 2015 with a 5km focal depth. The earthquake of ML1.4 occurred at 12:05 on December 15 with a 5km focal depth. The two earthquakes are located between the airgun source and Station 53262. The distance between the station and the airgun launching station is 4.5km. In this paper, the airgun signals recorded by stations close to the two earthquake epicenters are selected. According to the ray theory, signal rays received by Station 53262 are determined to have passed through the earthquake occurrence area, and those received by Station 53272 are not considered having passed through the source region.

The theoretical model experiments carried out above indicate that the cross-correlation delay method is more reliable than the dynamic programming method. Therefore, the cross-correlation delay method is usually used to calculate the wave velocity variation and analyze its possible correlation with the earthquake.

It can be seen from Fig. 5 that the actual active source signal is more complicated than the theoretical signal, and the phase is difficult to identify. The reason is probably that the process of triggering the seismic wave by the airgun source relies on the bubble oscillation process at the bottom of the water. The airgun wavelet, including pressure pulses and bubble pulses, is much more complicated than the Rake wavelet in the experiment (Tang Jie et al., 2009). In the process, three relatively stable seismic phases are selected as the object to calculate the wave velocity change. For Station 53262, the three phases have the arrival time periods of 2.5-3.0s, 3.0-3.5s and 3.5-4.2s, respectively. The seismic phase of Station 53272 is selected with reference to the corresponding seismic phase of Station 53262. The calculated delay result is shown in Fig. 5.

 Fig. 4 Locations of stations, airgun sources and earthquakes F1 is the Chenghai fault. F2 is the Honghe fault

 Fig. 5 Relative arrival time changes before and after the earthquake (a) is the template waveforms of Station 53262.(b), (c) and (d) are the wave-velocity changes of the three phases of Station 53262.(e), (f) and (g) are the wave-speed changes of the corresponding phases of Station 53272. The numbers in the figure indicate the sequence number of the events, wherein the adjacent events differed by 1 hour. Red and green colors indicate the relative changes in delay

First, according to the calculation in the theoretical model experiment, a rectangle which is 3km long and 1km wide is used, and the P-wave disturbing wave velocity in the body is 9km/s, the maximum delay that may occur in the ray path is 0.01s, much larger than the current measurement result. Usually, the change in travel time 10-5s can be produced even by a magnitude 1.0 earthquake (Zhou Qingyun et al., 2018), which is significantly less than the measurement results herein. The occurrence of this situation may be related to the excessive setting of parameters such as disturbance scale and medium change in the theoretical model. There are also other interference factors including solid tide.

Furthermore, there are obvious daily variations in wave velocity changes, which is consistent with previous research results.

Again, from the calculation results of 53262 units, it can be seen that the wave speed varies from falling to rising at 5 hours before the ML1.1 earthquake, which reaches its maximum 0.02% when the earthquake occurs. This phenomenon is more obvious before the earthquake. The first two phases show that there is a maximum time delay 5 hours before the earthquake, and the relative standard time delays are 9.50×10-4s and 8.07×10-4s respectively, while there is no obvious change with the third phase. This may be because the ray path of the third phase is not related to the area affected by the earthquake with ML1.4. Also, the adjacent Station 53272 does not show and drop before the earthquake with ML1.1. The three phases received at Station 53272 before the earthquake with ML1.4 showed the process of the wave speed falling to the rebound, but the change range was smaller than Station 53262. As shown in Fig. 4, the distance from Station 53262 to the airgun source is 4.5km, and the distance from Station 53272 to the airgun source is 5.2km. The two stations are 3.8km apart from each other. Besides, the two earthquakes are on the connecting line between the airgun source and Station 53262 and are farther away from Station 53272. The magnitudes of the earthquake differ by 0.3, which is the possible reason that the wave velocity of Station 53272 changed before the earthquake with ML1.4, but did not change before the earthquake with ML1.1.

Finally, when using the cross-correlation interferometry to measure the wave velocity, the main factors affecting the measurement results are the reservoir water level where the airgun source is located, the local temperature, the performance of the observation system, the calculation error, and the solid tide. In this experiment, because of its 170 hours duration time, the water level of the reservoir where the airgun is located was 22.67m at 08:00:00 on December 12, 2015, and 22.76m at 20:00:00 on December 18, 2015, during which the water level rose first and fell subsequently. The maximum value was 22.83m, and the change process was stable with a small amplitude. Zhou Qingyun (2018) calculated that the 3×10-4s time delay can be caused by the change of reservoir water level per 1m. During the experiment, the temperature in Binchuan area was stable, and the relationship between temperature and travel time was about 0.5μs/℃ (Niu Fenglin et al., 2008), much smaller than the current measurement results. If the performance of the observation system and the calculation error were not considered, the main interference factor in this experiment should be solid tide. According to Wang Weitao et al. (2017), the change in travel time has a great response to solid tide, which was confirmed by this experiment. The trend of wave velocity before the two earthquakes is not consistent with the trend of solid tides. To some extent, the measured wave velocity changes before the two earthquakes are reliable.

4 CONCLUSION

In this paper, focusing on the needs of the dynamic monitoring of underground media, large-capacity airgun source was taken as the research object, and the method of calculating the relative change of wave velocity using artificial source was analyzed.Consequently, several new insights were obtained:

(1) In theory, the dynamic programming method can tolerate a greater degree of deformation of the waveform than the interference method.

(2) In the actual calculation, the value of the wave velocity obtained by using the cross-correlation delay method to calculate the airgun source data is an order of magnitude larger than that of the dynamic gun method, but the wave velocity change trend is basically the same, which indicates the interference-based cross-correlation delay method for wave velocity measurement is more reliable.

(3) Under appropriate conditions, the Binchuan active source airgun signal is able to capture the process of small wave velocity reduction and recovery of an earthquake small as ML1.1.

REFERENCES
 Anderson K. R., Gaby J. E.Anderson K. R., Gaby J. E. Dynamic waveform matching[J]. Information Sciences, 1983, 31(3): 221-242. DOI:10.1016/0020-0255(83)90054-3 Chen Yong, Wang Baoshan, Ge Hongkui, Xu Ping, Zhang WeiChen Yong, Wang Baoshan, Ge Hongkui, Xu Ping, Zhang Wei. Proposed of transmitted seismic stations[J]. Advances in Earth Science, 2007, 22(5): 441-446 (in Chinese with English abstract). Cheng Xin, Niu Fenglin, Wang BaoshanCheng Xin, Niu Fenglin, Wang Baoshan. Coseismic velocity change in the rupture zone of the 2008 MW 7.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 Hale D.Hale D. Dynamic warping of seismic images[J]. Geophysics, 2013, 78(2): S105-S115. DOI:10.1190/geo2012-0327.1 Hall S.A.Hall S.A. A methodology for 7d warping and deformation monitoring using time-lapse seismic data[J]. Geophysics, 2006, 71(4): 021-031. Keysers D., Unger W.Keysers D., Unger W. Elastic image matching is NP-complete[J]. Pattern Recognition Letters, 2003, 24(1/3): 445-453. Li Le, Chen Qifu, Cheng Xin, Niu FenglinLi Le, Chen Qifu, Cheng Xin, Niu Fenglin. Spatial clustering and repeating of seismic events observed along the 1976 Tangshan fault, north China[J]. Geophysical Research Letters, 2007, 34(23): L23309. DOI:10.1029/2007GL031594 Lin Jianmin, Wang Baoshan, Ge Hongkui, Chen Qifu, Chen YongLin 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 Lanfeng, Wei Xiucheng, Huang Zhongyu, Wang Lu, Li LihuaLiu Lanfeng, Wei Xiucheng, Huang Zhongyu, Wang Lu, Li Lihua. PP-and PS-wave automatic time match based on dynamic image warp[J]. Oil Geophysical Prospecting, 2015, 50(4): 626-632 (in Chinese with English abstract). Liu Mian, Luo Gang, Wang HuiLiu Mian, Luo Gang, Wang Hui. The 2013 Lushan Earthquake in China tests hazard assessments[J]. Seismological Research Letters, 2014, 85(1): 40-43. DOI:10.1785/0220130117 Luo Guichun. Precise Measurement of Seismic Velocity and Velocity Variation by Correlated Detection[D]. Beijing: Institute of Earthquake Science, CEA, 2006 (in Chinese with English abstract). Nakata N., Snieder R.Nakata N., Snieder R. Near-surface weakening in Japan after the 2011 Tohoku-Oki earthquake[J]. Geophysical Research Letters, 2011, 38(17): L17302. 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 Pei Shunping, Niu Fenglin, Ben-Zion Y., Sun Quan, Liu Yanbing, Xue Xiaotian, Su Jinrong, Shao ZhigangPei Shunping, Niu Fenglin, Ben-Zion Y., Sun Quan, Liu Yanbing, Xue Xiaotian, Su Jinrong, Shao Zhigang. Seismic velocity reduction and accelerated recovery due to earthquakes on the Longmenshan fault[J]. Nature Geoscience, 2019, 12(5): 387-392. DOI:10.1038/s41561-019-0347-1 Pishchulin L. Matching Algorthms for Image Recognition[D]. Rheinisch-Westf alischen: Technischen HochSchule Aachen, 2010. Sakoe H., Chiba S.Sakoe H., Chiba S. Dynamic programming algorithm optimization for spoken word recognition[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1978, 26(1): 43-49. DOI:10.1109/TASSP.1978.1163055 Silver P. G., Daley T. M., Niu Fenglin, Majer E. L.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 Tang Jie, Wang Baoshan, Ge Hongkui, Chen YongTang Jie, Wang Baoshan, Ge Hongkui, Chen Yong. Study of experiment and simulation of large volume air-gun in deep structures exploration[J]. Earthquake Research in China, 2009, 25(1): 1-10 (in Chinese with English abstract). 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, Wang Weitao, Ge Hongkui, Xu Ping, Wang BinWang 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 Bin. The Experiment Study on Measurement of Seismic Velocity Variation Using Different Seismic Sources[D]. Hefei: University of Science and Technology of China, 2009 (in Chinese with English abstract). Wang Bin, Li Xiaobin, Liu Zifeng, Yang Jun, Yang Runhai, Wang Baoshan, Yang WeiWang Bin, Li Xiaobin, Liu Zifeng, Yang Jun, Yang Runhai, Wang Baoshan, Yang Wei. The source and observation system of Binchuan earthquake signal transmitting seismic station and its preliminary observation results[J]. Earthquake Research in China, 2016, 32(2): 193-201 (in Chinese with English abstract). Wang Weitao. Using Active Source to Study the Velocity Character of Media in Regional Scale[D]. Hefei: University of Science and Technology of China, 2009 (in Chinese with English abstract). Wang Weitao, Wang Baoshan, Ge Hongkui, Chen Yong, Yuan Songyong, Yang Wei, Li YijinWang 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). 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): 514-524 (in Chinese with English abstract). DOI:10.3969/j.issn.1000-0666.2017.04.002 Wang Zijie, Zhao Pengda. Primary results in studying geologic structural effect on earthquake ground motion with supersonic modelling technique[J].Geoscience——Journal of Graduate School, China University of Geosciences, 1997(3): 339, 353 (in Chinese with English abstract). Xu Z. J., Song XiaodongXu Z. J., Song Xiaodong. Temporal changes of surface wave velocity associated with major Sumatra earthquakes from ambient noise correlation[J]. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(34): 14207-14212. DOI:10.1073/pnas.0901164106 Yang Xiaosong, Ma JinYang Xiaosong, Ma Jin. Continental lithosphere decoupling:implication for block movement[J]. Earth Science Frontiers, 2003, 10(S1): 240-247 (in Chinese with English abstract). Yang Wei, Wang Baoshan, Ge Hongkui, Wang Weitao, Chen YongYang Wei, Wang Baoshan, Ge Hongkui, Wang Weitao, Chen Yong. 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 Wei. Finite Difference Seismic Wave Modelling in 3D Heterogeneous Media with Surface Topography and its Implementation in Strong Ground Motion Study[D]. Beijing: Peking University, 2006 (in Chinese with English abstract). Zhou Longquan, Liu Guiping, Ma Hongsheng, Hua WeiZhou Longquan, Liu Guiping, Ma Hongsheng, Hua Wei. Monitoring crustal media variation by using repeating earthquakes[J]. Earthquake, 2007, 27(3): 1-9 (in Chinese with English abstract). Zhou Qingyun, Liu Zifeng, He SugeZhou Qingyun, Liu Zifeng, He Suge. Influence factors of velocity variation of Airgun seismic waves and travel time changes related to earthquakes[J]. Earthquake, 2018, 38(3): 144-157 (in Chinese with English abstract).