Earthquake location is one of the most primary and classic problems, mainly focusing on improving the precision of earthquake location to the utmost. There are many factors affecting the precision of earthquake location such as the distribution of the record stations and the precision of the seismic phase picking and the accuracy of the speed model, for which errors about the speed model and seismic phase picking have a huge effect. Based on the double-difference method of waveform cross-correlation techniques, two errors could be reduced (Huang Yuan, 2008), with the improvement of the range of digital seismic waveform application. Using cross-correlation techniques to rectify the travel time of P-wave and S-wave could mitigate the influence of error in reading, and to a large degree, decrease errors caused by directly reading the seismic phase, solve the distributive problem of earthquake location caused by the inaccurate speed model and finally elevate the precision of earthquake location.
Cross-correlation analyses could make use of two methods, the time domain and frequency domain in calculation (Poupinet G. et al., 1984). The operation of the time domain method is simple and widely popular in application (Schaff D. P. et al., 2004, 2005; Waldhauser F. et al., 2000). Also, in terms of the limitation of cross-correlation calculation of the vertical component, multi-channel correlation inspection function of the time domain is developed and applied to calculate the travel time difference of waveform cross-correlation (Wang Qingdong et al., 2015).
The selection of the threshold of cross-correlation coefficient in the time domain needs some experience. If the threshold is set much higher, there is much less data available, so it is difficult to meet the requirements of high-precision location. Otherwise, some unreal data about seismic phase could be acquired. A relatively safer method is to use the bi-spectrum method to verify the cross-correlation coefficient results (Du Wenxuan et al., 2004). Through inspection of the cross-correlation coefficient, available information could be strengthened more effectively and cross-correlation data confidence could be elevated. The bi-spectrum method could repress Gaussian cross-correlation noise effectively, which could simultaneously calculate the two time delays of raw and band-pass-filtered waveforms in the third spectral domain, which could be used to examine the credibility of th ecross-correlation coefficient. With the characteristics of high precision, it has been put into use to some extent since it was proposed (Bannister S. et al., 2011; Bourguignon S. et al., 2015; Hansen S. E. et al., 2013; Zhang Guangwei et al., 2015; Zhao Cuiping, 2006).
The passage first uses the bi-spectrum software package (BCSEIS) to perform cross-correlation analysis of waveform data on seismic events in the Three Gorges Reservoir and precisely locate seismic events using double-difference location to compare locations by varied data to acquire the high-precision location of the Three Gorges Reservoir micro quake.
1 PROCESSING OF EARTHQUAKE WAVEFORM DATAThe Three Gorges earthquake observation network is composed of 26 earthquake observation stations, including 21 UK L-22 three-direction short-period speed pendulums and 5 American Guralp CMG-3ESPC wide-brand seismometers giving the DAQ of Reftek130 and a sampling rate of 200Hz (Ma Wentao et al., 2010). For the period from March, 2009 to December, 2010, there were a total of 5, 275 reports about seismic events and their waveform data.
In this research, we first conduct band-pass filtering of 1.0Hz-10Hz seismic waveform data, acquiring raw and band-pass-filtered waveform data, then conduct the cross-correlation analysis of the waveform cut out before and after the P and S waves arrive. We have 128 sampling points for the P-wave time window length, 30 sampling points before the arrival of the P-wave and 97 sampling points after its arrival. We have 192 sampling points for the S-wave time window length, 50 sampling points before the S-wave arrival and 141 sampling points after its arrival. Fig. 1 shows one event-pair waveform record, two seismic events numbered 242 and 244, in which the record from the original time of the earthquake to the end of S-wave is kept without other waveforms used to calculate the waveform cross-correlation coefficient. Detailed information about events 242 and 244 could be seen in Table 1. The calculation of waveform cross-correlation shows that event-pairs are highly similar, with two seismic event waveform cross-correlation coefficients higher than 0.9 and seismic phase time delays lower than 0.04s. Comparing the time difference of arrival of seismic phases in the seismic catalogue and the time delay difference of waveform cross-correlation, the waveform cross-correlation time delay could alleviate errors caused by seismic phase pickup errors and acquire more accurate data about the time difference of arrival, proving the waveform cross-correlation coefficient obtained from bi-spectrum verification is relatively precise.
In the process of bi-spectrum calculation, we use three cross-correlation coefficient thresholds to examine calculated results; central range coefficient CC^{lim(u)}, lower limit CC^{lim(l)} and upper limit CC^{lim(u)}. Usually, CC^{lim(l)} ≤CC^{lim} ≤CC^{lim(u)}, in which CC^{lim} is equal to the threshold in the calculation of traditional cross-correlation time delays. In practice, we make use of the following principles;
(1) If in the event pair (i, j), the maximum cross-correlation coefficient CC_{ij}^{max} ≥CC^{lim(u)}, we would accept cross-correlation values larger than from all stations, where cross-correlation coefficients are larger than CC^{lim(l)}. Otherwise, if they are smaller than CC^{lim(l)}, we can't accept them. In other words, regarding the seismic event pair with a highly similar waveform, we will choose the cross-correlation time delay results with cross-correlation coefficients larger than CC^{lim(l)} assessed by bi-spectrum verification. Even if we could acquire a cross-correlation coefficient smaller than CC, some researchers will not adopt the results smaller than the standard thresholds
(2) If CC^{lim} ≤CC _{ij}^{max} ≤CC^{lim(u)}, we only examine data from stations where cross-correlation coefficients are higher than CC^{lim}.
(3) If CC _{ij}^{max} is smaller than CC^{lim}, the cross-correlation time delays assessed results of seismic event pairs will be given up.
We use the following parameter settings to examine cross-correlation time delays: CC^{lim(l)} =0.5, CC^{lim} =0.7, CC^{lim(u)} =0.9, Δ^{lim} =10, among them, Δ^{lim} is the maximum lags between the band-pass filtering and the raw waveform data, using bi-spectrum verification. In order to compare with the traditional method of estimating cross-correlation coefficient thresholds, we select the cross-correlation coefficient threshold as 0.7, acquiring the traditional cross-correlation time delay data shown in Fig. 3. With the help of the traditional method of estimating cross-correlation coefficient thresholds, we get 1, 492, 711 P-wave seismic phase pairs and 749, 109 S-wave seismic phase pairs. With the help of bi-spectrum examination of cross-correlation time delays, there is a total of 1, 375, 372 P-wave seismic phase pairs and 749, 109 S-wave seismic phase pairs, among which there are 1, 273, 958 P-wave seismic phase pairs with cross-correlation coefficients larger than 0.7 and 596, 212 S-wave seismic phase pairs. Through bi-spectrum verification, the final seismic phase pairs are relatively smaller than the traditional cross-correlation seismic phase data. The main reason is that bi-spectrum verification eliminates cross-correlation data with relatively high time delay waveforms. The S-wave seismic phase pairs are less than the P-wave seismic phase pairs, the reason for which being that differentiating S waves is much harder than P waves. Also, S-waves are interfered with by P coda waves, which reduces the quality of the cross-correlation coefficient calculated within the time window of the S-wave.
In double difference location or double difference tomography, the maximum distance of seismic pairs affects the retrieved results greatly. The less the value is, the less the distance between seismic pairs establishing some relationship, and the more precise the localization is. However, the value is too small, leading to less seismic pairs establishing some relationship (Huang Yuan et al., 2006). Before double difference location or double difference tomography, it is necessary to assess the distance of seismic pairs. Normally, we can calculate the distance of seismic pairs and draw a depth-distribution histogram in order to decide the appropriate seismic distance as quickly as possible. We calculate the distance of seismic events within our study Fig. 4, which could indicate that it is rational to choose 20km as the distance of seismic pairs, consistent with the maximum distance of the strong connection given by the program ph2dt. However, whether the seismic distance matches the distance of cross-correlation seismic pairs requires further analysis. According to waveform cross-correlation data acquired by bi-spectrum verification, we collect event pairs "hypocentral distances and corresponding waveform cross-correlation coefficients from all stations and acquire the distributional graph of cross-correlation coefficients and event pairs" seismic source distances (Fig. 5). As a result, simply setting the distance threshold between seismic pairs as 20km will reduce the quality of the cross-correlation data. The localizing results are influenced by lateral heterogeneity, so it is necessary to reduce the threshold of the distance between the pairs, where we select 6km to ensure the consistency of most of the cross-correlation coefficients and distances, and to ensure the efficiency of cross-correlation data and the accuracy of localization.
Based on data and statistics from encrypted stations of the Institute of Geology, China Earthquake Administration between March, 2009 and December, 2010, there are a total of 5, 275 seismic events. Referring to the earthquake catalogue, they amount to 43, 538 items on P-wave travel time and 43, 385 S-wave items with 7 clear records for every seismic station. To exclude reading errors in seismic phases and other interfering information, we apply the least square fitting travel-time curve to winkle out relatively worse data and meanwhile examine the final results using the Wadati curve (namely the curve T_{S} -T_{P}) (Lin Yong et al., 2014). The Fig. 6 shows that there is a series of interfering information and then we process the data using the method of the least square fitting to eliminate relatively worse seismic phases, leading to relatively better information of seismic events. After election, we acquire 5, 244 seismic events, including data of seismic phases of wave P 41423 and wave S 42354. Fig. 7 shows the distribution of 5244 seismic events including information about longitudes, latitudes and depths. As a result, seismic events happen relatively intensively, mainly within the range of depth of over 8km. Seismic events are dispersed on the periphery of the Yangtze River system.
According to the results of the extant velocity model in the Three Gorges Reservoir Area (Li Qiang et al., 2009; Zhao Xu et al., 2007), the elected velocity models' horizontal depths are 0, 2, 5, 8, 11, 14 and 20km with corresponding speeds of wave P 4.8, 5.4, 5.65, 5.8, 6, 6.15 and 6.5 km/s. Based on fitting results, the value of v_{P} /v_{S} of the Fig. 6 is set as 1.73.
3 HIGH-PRECISION LOCALIZING RESULTS AND ANALYSES OF THE THREE GORGES RESERVOIR EARTHQUAKESAccording to the results of BCSEIS waveform cross correlation analysis and seismic catalogues, we use double difference localization software to test three groups of seismic events in the Three Gorges Reservoir area:the first group uses only the seismic catalogue to locate (CAT); the second group uses standard waveform cross-correlation methods to acquire waveform time delay and time difference of arrival of the seismic catalogue to locate (CC + CAT). The third group uses bi-spectrum verification to obtain waveform delay data and the seismic catalogue's time difference of arrival to locate (WCC + CAT). In order to compare the differences in localization, the same parameters and velocity models were set for the three sets of test data, to finally gain the seismic planar position (Fig. 8) and the depth error distribution of different data types (Fig. 9). Fig. 8 shows that the localizations of the second and third groups are significantly better than those of the first set of data, with smaller errors and predominant depths (Fig. 9).
The epicenters after relocation are more concentrated and linearly distributed, showing the trend of inward convergence of earthquake events. There are three obvious linear distributions almost from east to west on both sides of Shenlongxi river in Badong with average depth of about 5km. The epicentral distribution is basically consistent with the stratigraphic trend, which further confirms the reservoir water storage from the Shenlongxi river leading to other underground rivers' infiltration and inducing the earthquakes (Ma Wentao et al., 2010). A M_{S} 5.1 earthquake occurred east of the seismic distributive line on December 16, 2013, and the microseismic localization from March, 2009 to December, 2010 indicates that the earthquake under the effect of reservoir water caused a larger earthquake activity that triggered a small-scale rupture through the east-west. In the western part of the Xietan area, there is a steep east-west microseismic activity zone, which may be related to the infiltration of the Yangtze River system.In the vicinity of the Xiangxi river in Zigui County, the earthquakes are distributed on the fault zone of the Xiannvshan mountain, and the depth of the source is extended down to 10km, revealing the maximum depth of reservoir water infiltration. This is where two M_{S} 4.0 earthquakes occurred, which may have also been caused by reservoir water infiltration arousing regional activities in the fault zone.
The localizing errors' statistics of the three different types of data (Fig. 10) show that the hypoDD positioning accuracy of the waveform cross-correlation delay using bi-spectrum verification and the time difference of the seismic catalogue are significantly higher than that of the other two kinds of data. The average error of hypocenters' localities is 3.2m in the east-west direction, 3.9m in the north-south direction and 6.2m in the vertical direction.
In this paper, the bi-spectrum method is used to analyze the seismic observation data of the Three Gorges Reservoir area from March, 2009 to December, 2010 in terms of waveform cross-correlation analysis. The high-quality bi-spectrum verified waveform cross-correlation delays data are obtained and the double difference method is used to relocate earthquakes to obtain high-precision positioning results.
Based on the third-order spectral domain, the bi-spectral waveform cross-correlation algorithm can derive higher quality waveform delay data than that by using the traditional threshold, and the distance of the cross-correlation seismic event pairs is more realistic, which can provide higher quality data for double difference or tomography.
According to the traditional threshold judgment the traditional cross-correlation data obtained by bi-spectral verification, double difference localizing programs are respectively used. The results show that the accuracy of localization of the waveform cross-correlation data is higher than that of other data. As a result, the average error of the east-west hypocentral location is 3.2m, 3.9m in the south-north direction and the vertical direction is 6.2m.
The results of the Three Gorges Reservoir area's seismic relocation show that the microseismic distribution along two sides of Shenlongxi river in the Badong reservoir area indicates three linear distributions in the east-west direction, which is consistent with the trends of the small-scale faults and carbonate strata in the neotectonic period. The reservoir water travels mainly along the cave or underground river to infiltrate and induce seismic activities, among which the seismic distributional line endpoint is the position of epicenter of the occurrence of the December 16, 2013 Badong, Hubei M_{S} 5.1 earthquake, indicating that there was an event of larger rock mass dislocation of a small-scale fault in the east-west direction under the effect of the reservoir water. Xietan's western east-west vertical stripe-shaped distribution shows the relationship between the seismic activities and the reservoir water infiltration in the Yangtze River. In the vicinity of the Xiangxi river in Zigui County, earthquakes are distributed on the fault zone of the Xiannvshan, and the depth of the source is extended down to 10km, which reveals the maximum depth of reservoir water infiltration. d
ACKNOWLEDGEMENT
Sincere thanks extended to reviewing experts who offered critical advice and GMT, the software applied to create most figures in this article (Wessel P. et al., 1995).
This paper was published in Chinese in the journal of Technology for Earthquake Disaster Prevention, Volume 12, Number 1, 2017.
Bannister S., Fry B., Reyners M., Ristau J., Zhang Haijiang. Fine-scale relocation of aftershocks of the 22 February M_{W}6.2 Christchurch earthquake using double-difference tomography[J]. Seismological Research Letters, 2011, 82(6): 839–845. DOI:10.1785/gssrl.82.6.839. |
Bourguignon S., Bannister S., Henderson C.M., Townend J., Zhang Haijiang. Structural heterogeneity of the midcrust adjacent to the central Alpine Fault, New Zealand:Inferences from seismic tomography and seismicity between Harihari and Ross[J]. Geochemistry, Geophysics, Geosystems, 2015, 16(4): 1017–1043. DOI:10.1002/2014GC005702. |
Du Wenxuan, Thurber C.H., Eberhart-Phillips D. Earthquake relocation using cross-correlation time delay estimates verified with the bi-spectrum method[J]. Bulletin of the Seismological Society of America, 2004, 94(3): 856–866. DOI:10.1785/0120030084. |
Hansen S.E., DeShon H.R., Moore-Driskell M.M., Al-Amri A.M.S. Investigating the P wave velocity structure beneath Harrat Lunayyir, northwestern Saudi Arabia, using double-difference tomography and earthquakes from the 2009 seismic swarm[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(9): 4814–4826. DOI:10.1002/jgrb.50286. |
Huang Yuan, Yang Jiansi, Zhang Tianzhong. Relocation of the Bachu-Jiashi, Xinjiang earthquake sequence in 2003 using the double-difference location algorithm[J]. Chinese Journal of Geophysics, 2006, 49(1): 162–169. |
Huang Yuan. Study on the application and development of the DD algorithm with cross correlation of waveform data in the earthquake location[J]. Recent Developments in World Seismology, 2008(4): 29–34. |
Li Qiang, Zhao Xu, Cai Jin'an, Liu Ruifeng, Long Guihua, An Yanru. P wave velocity structure of upper and middle crust beneath the Three Gorges reservoir dam and adjacent regions[J]. Science in China (Series D:Earth Sciences), 2009, 52(4): 567–578. DOI:10.1007/s11430-009-0047-6. |
Lin Yong, Ma Wentao. Matlab-based seismic data processing and analysis of Three Gorges Reservoir[J]. Technology for Earthquake Disaster Prevention, 2014, 9(3): 447–453. |
Ma Wentao, Xu Changpeng, Li Haiou, Yuan Jingli, Xu Xiwei, Zhang Xindong, Zhang Lanfeng. Intensive observation of reservoir-induced earthquake and preliminary analysis on the causes of earthquakes in Three Gorges Reservoir[J]. Seismology and Geology, 2010, 32(4): 552–563. |
Poupinet G., Ellsworth W.L., Frechet J. Monitoring velocity variations in the crust using earthquake doublets:An application to the Calaveras Fault, California[J]. Journal of Geophysical Research:Solid Earth, 1984, 89(B7): 5719–5731. DOI:10.1029/JB089iB07p05719. |
Schaff D.P., Bokelmann G.H.R., Ellsworth W.L., Zanzerkia E., Waldhauser F., Beroza G.C. Optimizing correlation techniques for improved earthquake location[J]. Bulletin of the Seismological Society of America, 2004, 94(2): 705–721. DOI:10.1785/0120020238. |
Schaff D.P., Waldhauser F. Waveform cross-correlation-based differential travel-time measurements at the Northern California Seismic Network[J]. Bulletin of the Seismological Society of America, 2005, 95(6): 2446–2461. DOI:10.1785/0120040221. |
Waldhauser F., Ellsworth W.L. A double-difference earthquake location algorithm:Method and application to the northern Hayward fault, California[J]. Bulletin of the Seismological Society of America, 2000, 90(6): 1353–1368. DOI:10.1785/0120000006. |
Wang Qingdong, Zhu Liangbao, Su Youjin, Wang Guangming. Double-difference relocation of the 7 September 2012 Yiliang earthquake and its aftershock sequence[J]. Chinese Journal of Geophysics, 2015, 58(9): 3205–3221. |
Wessel P., Smith W.H.F. New version of the generic mapping tools[J]. Eos, Transactions American Geophysical Union, 1995, 76(33): 329. |
Zhang Guangwei, Lei Jianshe. Mechanism of the 2011 Tengchong, Yunnan, M_{S}5.2 double earthquakes[J]. Chinese Journal of Geophysics, 2015, 58(4): 1194–1204. |
Zhao Cuiping. Seismological Dtudies on the Characteristics of Jiashi Source Region from 1977 to 2003[D]. Doctoral thesis. Beijing: Institute of Geophysics, China Earthquake Administration, 2006 (in Chinese with English abstract). |
Zhao Xu, Li Qiang, Cai Jin'an. On minimum 1D velocity model applied in Three Gorges Reservoir area[J]. Journal of Geodesy and Geodynamics, 2007, 27(S1): 1–7. |