2. University of Science and Technology of China, Hefei 230000, China
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 earthquaketriggered 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 noiserelated 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 highprecision data of wave velocity was obtained from rock physics experiments, but the measurement of highprecision 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 highprecision observation systems, airgun has received much attention as a important source so that largecapacity 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 crosscorrelation 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 CrossCorrelation Delay MethodDue 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 crosscorrelationbased 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 subsampling 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 powerlimited 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 ProgrammingThe 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 highprecision 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(k1) \leqslant i(k) & \text { and } j(k1) \leqslant j(k) \\ \end{aligned} $  (10) 
$ \begin{aligned} i(k)i(k1) & \leqslant 1 \text { and } j(k)j(k1) \leqslant 1 \end{aligned} $  (11) 
that is
$ c(k1)=\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 ModelThere 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 twodimensional 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 mediumdisturbing 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.
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 Swave is greater than that on Pwave (Fig. 2).
2.2 Analysis of Experimental ResultsThe crosscorrelation 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, timelapse 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 crosscorrelation 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.
As shown in Fig. 3, the crosscorrelation 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 crosscorrelation 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 crosscorrelation method based on interference are too small. In addition, the calculation results of the Pwave 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 crosscorrelation delay method, the dynamic programming method has a more obvious amplification effect and is more sensitive to small disturbances. Correspondingly, the crosscorrelation method has strong stability.
(3) The wave calculation results obtained by the two methods show that the change trend of Pwave calculation is almost the same and the calculation result difference of the Swave 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 Swave 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 crosscorrelation 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 AREADuring December 1216, 2015, the Binchuan Seismic Signal Transmitting System has continuously conducted a 4day 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 M_{L}1.1 occurred at 20:46 on December 12, 2015 with a 5km focal depth. The earthquake of M_{L}1.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 crosscorrelation delay method is more reliable than the dynamic programming method. Therefore, the crosscorrelation 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.53.0s, 3.03.5s and 3.54.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.
First, according to the calculation in the theoretical model experiment, a rectangle which is 3km long and 1km wide is used, and the Pwave 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^{5}s 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 M_{L}1.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^{4}s and 8.07×10^{4}s 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 M_{L}1.4. Also, the adjacent Station 53272 does not show and drop before the earthquake with M_{L}1.1. The three phases received at Station 53272 before the earthquake with M_{L}1.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 M_{L}1.4, but did not change before the earthquake with M_{L}1.1.
Finally, when using the crosscorrelation 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^{4}s 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 CONCLUSIONIn this paper, focusing on the needs of the dynamic monitoring of underground media, largecapacity 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 crosscorrelation 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 interferencebased crosscorrelation 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 M_{L}1.1.
Anderson K. R., Gaby J. E.Anderson K. R., Gaby J. E. Dynamic waveform matching[J]. Information Sciences, 1983, 31(3): 221242. DOI:10.1016/00200255(83)900543 
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): 441446 (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 M_{W} 7.9 Wenchuan earthquake observed from ambient seismic noise[J]. Bulletin of the Seismological Society of America, 2010, 100(5B): 25392550. DOI:10.1785/0120090329 
Hale D.Hale D. Dynamic warping of seismic images[J]. Geophysics, 2013, 78(2): S105S115. DOI:10.1190/geo20120327.1 
Hall S.A.Hall S.A. A methodology for 7d warping and deformation monitoring using timelapse seismic data[J]. Geophysics, 2006, 71(4): 021031. 
Keysers D., Unger W.Keysers D., Unger W. Elastic image matching is NPcomplete[J]. Pattern Recognition Letters, 2003, 24(1/3): 445453. 
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): 19 (in Chinese with English abstract). 
Liu Lanfeng, Wei Xiucheng, Huang Zhongyu, Wang Lu, Li LihuaLiu Lanfeng, Wei Xiucheng, Huang Zhongyu, Wang Lu, Li Lihua. PPand PSwave automatic time match based on dynamic image warp[J]. Oil Geophysical Prospecting, 2015, 50(4): 626632 (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): 4043. 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. Nearsurface weakening in Japan after the 2011 TohokuOki 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): 204208. DOI:10.1038/nature07111 
Pei Shunping, Niu Fenglin, BenZion Y., Sun Quan, Liu Yanbing, Xue Xiaotian, Su Jinrong, Shao ZhigangPei Shunping, Niu Fenglin, BenZion 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): 387392. DOI:10.1038/s4156101903471 
Pishchulin L. Matching Algorthms for Image Recognition[D]. RheinischWestf 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): 4349. 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 crosswell seismic travel time for stressinduced changes[J]. Bulletin of the Seismological Society of America, 2007, 97(1B): 281293. 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 airgun in deep structures exploration[J]. Earthquake Research in China, 2009, 25(1): 110 (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): 249256 (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): 193201 (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): 223233 (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): 514524 (in Chinese with English abstract). DOI:10.3969/j.issn.10000666.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): 1420714212. 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): 240247 (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): 399410 (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): 19 (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): 144157 (in Chinese with English abstract). 