Earthquake Reaearch in China  2017, Vol. 31 Issue (3): 368-380
Application of the Double-difference Relocation Method Combined with Waveform Cross-correlation on the Three Gorges Reservoir Seismicity
Luo Jiahong, Ma Wentao, Li Chunzheng     
Key Laboratory of Active Tectonics and Volcano, Institute of Geology, CEA, Beijing 100029, China
Abstract: In this paper, we use the double difference location method based on waveform cross-correlation algorithm for precise positioning of the Three Gorges Reservoir (TGR) earthquakes and analysis of seismic activity. First, we use the bi-spectrum cross-correlation method to analyze the seismic waveform data of TGR encrypted networks from March, 2009 to December, 2010, and evaluate the quality of waveform cross-correlation analysis. Combined with the waveform cross-correlation of data obtained, we use the double difference method to relocate the earthquake position. The results show that location precision using bi-spectrum verified waveform cross-correlation data is higher than that by using other types of data, and the mean 2 sig-error in EW, NS and UD are 3.2m, 3.9m and 6.2m, respectively. For the relocation of the Three Gorges Reservoir earthquakes, the results show that the micro-earthquakes along the Shenlongxi river in the Badong reservoir area obviously show the characteristics of three linear zones with nearly east-west direction, which is in accordance with the small faults and carbonate strata line of the neo-tectonic period, revealing the reservoir water main along the underground rivers or caves permeated and induced seismic activity.The stronger earthquakes may have resulted from small earthquakes through the active layers.
Key words: Bi-spectrum verification     Waveform cross-correlation     Double-difference location     The Three Gorges Reservoir    

INTRODUCTION

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 DATA

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

Fig. 1 Waveforms of event pairs from the Three Gorges Reservoir (TGR) intensive seismic observation station (The event pairs in Table 1) (a) The waveform record of Earthquake 242.(b) The waveform record of Earthquake 244

Table 1 Detailed information of event pairs shown in Fig. 1

In the process of bi-spectrum calculation, we use three cross-correlation coefficient thresholds to examine calculated results; central range coefficient CClim(u), lower limit CClim(l) and upper limit CClim(u). Usually, CClim(l)CClimCClim(u), in which CClim 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 CCijmaxCClim(u), we would accept cross-correlation values larger than from all stations, where cross-correlation coefficients are larger than CClim(l). Otherwise, if they are smaller than CClim(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 CClim(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

Fig. 2 Waveform analysis results of event pair 1 on station DTP (a) and LPT (b) The black waveform represents waveform 242, the red waveform represents waveform 244, P in the first two lines represents travel-time of P-wave of the earthquake catalogue of the station; CC in the third line represents waveform cross-correlation coefficient; tdt as the time-difference of arrival of waveform P of waveform cross-correlation, tdiff as time-difference of waveform P of the earthquake catalog

(2) If CClimCC ijmaxCClim(u), we only examine data from stations where cross-correlation coefficients are higher than CClim.

(3) If CC ijmax is smaller than CClim, 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: CClim(l) =0.5, CClim =0.7, CClim(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.

Fig. 3 Histogram of two kind of cross-correlation coefficients in the Three Gorges Reservoir Area (a) The statistical chart of traditional cross-correlation analysis of P waves; the values in the bracket in Y as the number; the statistical interval of cross-correlation coefficients is 0. 01. (b) The statistical chart of cross-correlation of P-waves using bi-spectrum verification; the first value in the bracket of Y is the number of cross-correlation coefficients larger than 0.7; the second is the total number; the three vertical lines indicate CClim(l) =0.5, CClim =0.7, CClim(u) =0.9.(c) The statistics of traditional cross-correlation analysis of S waves, the rest equal to (a). (d) The statistical chart of cross-correlation analysis of S-wave using bi-spectrum verification, the rest equal (b)

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.

Fig. 4 Histograms of event pair hypocentral separation distance for the event catalog

Fig. 5 Distribution of event pairs hypocentral separation distance and cross-correlation coefficient for all stations (a) The distribution of P waves data. (b) The distribution of P-wave cross-correlation coefficients and sources' distances density; X represents sources' distances of same seismic pairs from all stations, Y as corresponding cross-correlation coefficients; the interval source distance as 2km, cross-correlation coefficient as 0.1 in the statistical chart of density; (c) and (d) as distributions of S waves, the rest equal to (a) and (b)
2 REALIZATION OF LOCALIZATION OF DOUBLE DIFFERENCE IN THE THREE GORGES RESERVOIR AREA

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 TS -TP) (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.

Fig. 6 The curves of traveltime and Wadati plot (a) The travel-time curve after using the least square fitting to eliminate relatively worse seismic phases. (b) The primitive curve of travel-time of seismic phases. (c) The TS -TP after eliminating relatively worse seismic phases, the red line represents fitting results of the straight line. (d) The primitive curve TS -TP

Fig. 7 The epicenter distribution map, depth profile and statistical chart of 5, 244 seismic events (a) The primitive hypocentral distribution of seismic events. (b)The cross-section map of hypocentral depths along the latitude. (c) The cross-section map of hypocentral depths along the longitude. (d) The statistical graph of hypocentral depths

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 vP /vS of the Fig. 6 is set as 1.73.

3 HIGH-PRECISION LOCALIZING RESULTS AND ANALYSES OF THE THREE GORGES RESERVOIR EARTHQUAKES

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

Fig. 8 Map of relocation results from different type data by HypoDD (a) The primitive distribution of hypocenters of seismic events. (b) The localizations of the seismic catalogue's HypoDD. (c)The localizations of the seismic catalogue's HypoDD combined with traditional cross-correlation analysis. (d) the localizations of the seismic catalogue's HypoDD combined with bi-spectrum cross-correlation analysis, the left upper number representing the number of relocating seismic events using different data,

Fig. 9 The distribution of depth and RMS from different data by HypoDD relocation results (a)Data only referring to the earthquake catalog. (b)Data referring to earthquake catalog and traditional cross-correlation analysis. (c)Cross-correlation analysis using bi-spectrum and localizations of the earthquake catalog HypoDD; the upper as the depth-distribution in longitudes; the lower as the error-distribution in longitudes

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

Fig. 10 The error statistics of different type data by HypoDD relocation (a) -(c) the localization of the earthquake catalog. (e) -(g) the localization of traditional cross-correlation analysis and the earthquake catalog. (h) -(j) the localization of waveform cross-correlation using bi-spectrum verification and the earthquake catalog
4 CONCLUSION

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

REFERENCES
Bannister S., Fry B., Reyners M., Ristau J., Zhang Haijiang. Fine-scale relocation of aftershocks of the 22 February MW6.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, MS5.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.