Earthquake Research in China  2019, Vol. 33 Issue (4): 617-631     DOI: 10.19743/j.cnki.0891-4176.201904011
ETAS Model Analysis on the Chang Island Earthquake Swarm in Shandong Province, China
WANG Peng1,2, WANG Baoshan1,3     
1. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. Shandong Earthquake Agency, Jinan 250102, China;
3. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China
Abstract: Influenced by the layout of seismic network and the location of earthquakes, earthquake catalogs are often incomplete; such incompleteness of earthquake catalogue directly affects the analysis of sequence activity characteristics. In this paper, the GPU-acceleration-based g template matching method is used to scan the continuous waveforms of Chang Island earthquake swarm in Shandong Province from February 9 to August 20, 2017. In total, 15, 286 earthquakes events were detected, which was more than 6 times compared with those in network catalogue and thus reduced the magnitude of completeness from 1.0 to 0.5. Based on the intergrated catalogue of earthquakes, the characteristics of Chang Island earthquake swarm were then analyzed using the Epidemic Type Aftershock Sequences (ETAS) model. The stochastic components in the ETAS model are used as a proxy for possible earthquake triggered by external forces (fluids). The results show that the proportion of earthquakes triggered by external forces of Chang Island swarm increases gradually (from 31.9% to 63.5%) and then decreases. The latter stage of swarm development is mainly affected by the self-excitation of earthquakes, suggesting that the fluids play an important role in the development of the Chang Island swarm. However, the triggering intensity of fluids to microseismicity is divergent in different periods, which may be related to the process of fluid permeation.
Key words: ETAS model     Chang Island swarm     Template matching     Magnitude of completeness     b-value     Fluid triggering    

INTRODUCTION

Since February 14, 2017, seismic activities have occurred continually in Chang Island area, Shandong Province. Until August, 2017, more than 2, 300 earthquakes have been recorded in the local catalog, including 43 ML≥3.0 events. The largest earthquake was the ML4.5 earthquake occurred on March 3, 2017 (Fig. 1). The frequency and intensity of the Chang Island earthquake swarm exceed those before the 1976 Tangshan MS7.8 earthquake in the same area, which is unwonted in the history. Chang Island area has always been considered as an information window of the seismicity investigation in North China (Su Luansheng, 1993). Therefore, detailed analysis of the characteristics of Chang Island earthquake swarm is conducive to understand the seismogenic mechanism of earthquakesand further assist estimate the seismic risk in this area.

Fig. 1 The location of the Chang Island swarm (red dots) and local seismic stations (triangles) Blue triangles are the stations used for earthquake detection. The red box in the inset indicates the study area of the main map

The Chang Island swarm is located at the eastern end of the NWW-oriented Zhangjiakou-Penglai fault zone. The faults in this area are complex and mainly consist of the NWW-oriented Penglai-Weihai fault and several NE-trending small faults (Fig. 1). Due to the difficulty of tectonic stress measurement as well as the complexity of the seismic triggering process, it is difficult to establish a physical model for accurate description. However, the essential rules can be explored from complex phenomena using statistical methods. Since Ogata (Ogata Y., 1988) proposed the Epidemic Type Aftershock Sequence (ETAS) model based on the "Omori-Utsu" law (Omori F., 1894; Utsu T., 1961), the model has received extensive attention and gained considerable development. With the support of ETAS model, the parameters of seismic sequence can be accurately estimated, the variation law of seismic activity can be effectively found (Ogata Y., 1992, 2001; Jiang Haikun et al., 2007; Jiang Changsheng et al., 2013(a), 2013(c), 2015, 2017), the relative quietness of seismic activity can be monitored (Zhuang Jiancang, 2000; Jiang Changsheng et al., 2013(b)), slight stress change can be detected (Helmstetter A. et al., 2003; Lei Xinglin et al., 2008(a), 2017; Jia Ke et al., 2014, 2018), and induced earthquakes related to fluid action can be identified (Hainzl S. et al., 2005; Lei Xinglin et al., 2008(b); Long Feng et al., 2010; Jiang Haikun et al., 2011, 2012; Peng Yajun et al., 2012)

The statistical seismic modelling method depends largely on the earthquake catalogue. The seismic sequence activity characteristics and the aftershock probability prediction of large earthquakes are usually affected by the magnitude of completeness (Mc) (Jiang Changsheng et al., 2013(c); Zhuang Jiancang et al., 2017). Also, the artificial setting of the magnitude of completeness may affect the estimation of ETAS model parameters (Jiang Haikun et al., 2012; Li Zhichao et al., 2014), which mainly depends on the monitoring ability of regional seismic network. In areas like Chang Island, the signal-to-noise ratios of seismic waveforms are usually low, and seismic stations are sparse or the distribution is not ideal. Thus, the clustering occurrence of earthquake sequences in a short time usually results in the omission of earthquake in the catalog. However, the reliability of analysis can be improved by using a more complete earthquake catalogue (Schaff D. P., 2008).

The template matching filtering technique based on waveform cross-correlation is suitable for microseismic detection (Gibbons S. J. et al., 2006), which has developed rapidly and has been widely used recently. It was applied for detecting tremors (Shelly D. R. et al., 2007), studying the aftershock sequence changes (Peng Zhigang et al., 2009; Hou Jinxin et al., 2017; Wang Peng et al., 2017), determining the seismogenic structure (Yang Hongfeng et al., 2009; Tan Yipei et al., 2016), analyzing the stress changes in the source region (Meng Xiaofeng et al., 2012), identifying repeating earthquakes (Ma Tengfei, 2015), and monitoring nuclear explosions (Zhang Miao et al., 2015).

In this paper, a GPU-based template matching method is used to detect the missing earthquakes of the Chang Island earthquake swarm. The statistical parameters of the ETAS model of the earthquake swarm are analyzed based on the completed earthquake catalogue, and the possible seismogenic mechanism of the Chang Island earthquake swarm is then discussed.

1 GEOLOGICAL BACKGROUND AND SEISMIC ACTIVITY IN THE STUDY AREA

The Chang Island earthquake swarm is located on the Penglai-Weihai fault zone, which is the southeast segment of the NWW-trending Zhangjiakou-Bohai fault zone (Xu Jie et al., 1998). It is composed of a series of NW-trending faults developed between Penglai and Weihai sea area in the northern Shandong Peninsula. The major part is located between the Chang Island and the Dazhu Island; most faults are normal and strike-slip faults, some of which are thrusting (Zheng Jianchang et al., 2018). The Penglai-Weihai fault zone has hosted many historical earthquakes. Among them, the Chang Island MS7.0 earthquake in 1948 occurs in the west of the fault zone and the Weihai MS6.0 in 1948 occurs in the east of the fault zone (Wang Zhicai et al., 2006).

From February 14 to August 30, 2017, an earthquake swarm with 2, 377 earthquakes (from the Shandong Seismic Network) has occurred in the Chang Island, including 1, 459 ML1.0-1.9 earthquakes, 271 ML2.0-2.9 earthquakes, 39 ML3.0-3.9 earthquakes and 43 ML≥4.0 earthquakes. Among them, the largest ML4.5 earthquake occurred on March 3. There are 16 seismic stations within 150km from the swarm, all are three-component broadband seismographs with a sampling rate of 100Hz. They are mainly distributed in the land area to south of the swarm. In the following, we select Chang Island Station (CHD), Beihuangcheng Station (BHC), Longkou Station (LOK), and Yantai Station (YTA), which are close to and well distributed to the earthquake swarm, to detect missing earthquakes by using the three-component continuous waveforms.

2 METHODS 2.1 Template Matching Method

The procedure of the template matching method used in this paper follows the work by Peng and Zhao (Peng Zhigang and Zhao Peng, 2009), and uses the similar processing steps as Hou and Wang (Hou Jinxin and Wang Baoshan, 2017). The procedure is shown in Fig. 2, which is roughly divided into four steps and is briefly described as follows:

Fig. 2 Template matching method calculation process

(1) Data preprocessing and template selection. Using the three-component records of the four stations mentioned above (blue triangles in Fig. 1), the mean values of the continuous waveform data from the period February 9 to August 30, 2017 are removed. Then, the data are detrended and band-pass filtered to 2-10Hz. In total, there are 302 earthquakes with ML≥2.0 and SNR>5 from Chang Island swarm selected as templates, and the P- and S-wave arrival times of the template earthquakes are obtained from the local catalog.

(2) The acquisition of cross-correlation coefficients. The template waveforms are then cross-correlated to the corresponding continuous waveforms. The correlating time window is 1s before and 5s after the P- and S-wave arrival times of the template. The cross-correlation coefficient of the four stations are then shifted and averaged.

(3) The detection of missing earthquakes.We calculate the Median Absolute Deviation (MAD) of averagedcross-correlation coefficients, and select 12 times of the MAD as threshold. When the cross-correlation coefficient is higher than the threshold, it is marked as a tentative event; the maximum value of the average correlation coefficient within 5s is set as a detected event.

(4) The determination of earthquake time and magnitude. Based on the assumption that the missing event has similar travel times as the corresponding template event, the origin time of the missing event is obtained. The ratio of maximum amplitude of the horizontal component of the missing event to the corresponding template event is used to estimate the magnitude of the missing earthquake.

2.2 ETAS Model

The conditional intensity function λ(t) of the ETAS model is (Ogata Y., 1988; Jiang Haikun et al., 2012):

$ \lambda \left( t \right) = \mu \left( t \right) + K\sum\limits_{{t_i} < t} {\frac{{{e^{\alpha \left( {{M_i} - {M_z}} \right)}}}}{{{{\left( {t - {t_i} + c} \right)}^p}}}} , $ (1)

The first item on the right μ(t) is the theoretical number of earthquakes per unit time after eliminating the clustering of earthquakes. It also indicates that the intensity of microseismic activity triggered by external force (fluid), which has no correlation with the self-excitation process between earthquakes (Hainzl S. et al., 2005). The second item describes the contribution of the ith earthquake to the occurrence of following earthquakes. i traverse all earthquakes in the earthquake sequence. Mz is the reference magnitude no less than the magnitude of completeness Mc. p represents the attenuation factor of seismic sequence: large p means fast attenuation, while small pmeans slow attenuation. In the mainshock-aftershock type earthquake sequence, when Mz takes the magnitude of the mainshock, the smaller α value has higher secondary aftershock excitation ability (Song Jin et al., 2009); when Mz takes the complete magnitude of Mc, α indicates the ability of triggering secondary aftershocks, and greater α indicates the stronger ability (Jiang Haikun et al., 2012). The initial physical implication of c is a very small positive number, ensuring that the denominator of the second term of equation (1) is not zero, the exact value of c is determined by the completeness of the data in a short time after the mainshock (Jiang Haikun et al., 2007) and the stress state of the source region (Narteau C. et al, 2009). k is the level of seismic activity associated with the magnitude of the mainshock and the b, p, α, c, μ parameters, where b is the scale factor in the G-R law (Song Jin et al., 2009).

One observed seismic sequence can be fitted with multiple statistical models mathematically, andmultiple sets of qualified parameters can be obtained. Ogata (Ogata Y., 1988) uses the maximum likelihood method to estimate the parameters of the ETAS model and uses the Akaike information criterion (AIC) to minimize the optimal model parameters (Akaike H., 1974). The AIC criterion not only considers the pros and cons of the model fitting of the observed data, but also punishes the behavior of increasing the model parameters in an unrestricted manner for improving the fitting degree. The log likelihood of the ETAS model are:

$ \lg L = \sum\limits_{i = 1}^{{N_{{\rm{AIC}}}}} {\lg \lambda \left( t \right) - \int\limits_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\lambda \left( t \right){\rm{d}}t} } , $ (2)
$ {\rm AIC} = - \lg L + 2k, $ (3)

In equation (2), tb and te are the beginning and end time used for AIC calculation, and NAIC is the number of earthquakes in the calculation time period. It is necessary to ensure that the earthquakes in the calculation period are basically complete. In formula (3), k is the number of parameters to be fitted. In this paper, k=5, and the model with minimum AIC is taken as the final model.

It can be seen from equation (1) that the seismic rate represented by ETAS model consists of two parts, the first term is independent of the aftershock excitation triggered by the external forces, while the second item is the seismic cluster caused by the "aftershock excitation aftershock". Therefore, through the separation of the two parts of the ETAS model, the ratio of the external-force-triggered earthquake Rb can be expressed by formula (4), and the proportion of self-excited earthquake is Rt=1-Rb, The relationship between the triggering of external factors (fluid action) and seismic activity can be investigated according to their respective proportions. (Jiang Haikun et al., 2011).

$ {{R}_b} = \frac{{\int_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\mu \left( t \right){\rm{d}}t} }}{{\int_{{t_{\rm{b}}}}^{{t_{\rm{e}}}} {\lambda \left( t \right){\rm{d}}t} }}, $ (4)
3 RESULTS 3.1 Missing Earthquake Detection Results

Using the template matching technique, we detected 15, 286 earthquakes from February 9 to August 20, 2017, including 302 self-detected templates and 14, 984 missing events. The seismic catalogue obtained by template matching method is hereinafter referred to as detection catalog. The event number in the detection catalog is more than 6 times compared with the number in Shandong Seismological Network Catalog. Fig. 3 is an example showing the detection of a missing ML2.6 earthquake. The seismic rates of the detection and network catalogues are shown in Fig. 4, from which we suggest that there are few earthquakes when magnitude ML > 2.5. The difference between the detection and network catalogues is related to the monitoring ability of the network in this area. The monitoring ability of one network can be described with the magnitude of completeness. To evaluate the magnitude of completeness, we used the maximum curvature method, and the corresponding magnitude with the maximum slope in the cumulative frequency-magnitude curve is usually selected as the magnitude of completeness(Mc)(Wiemer S. et al., 2000). In practical applications, this magnitude usually corresponds to the magnitude of the non-cumulative frequency-magnitude distribution with the largest seismic frequency (Huang Yilei et al., 2016). The estimated completeness magnitudes of the detection catalogue and the network catalogue are 0.5 and 1.0, respectively. The cumulative frequency-magnitude curve is then fitted using the maximum likelihood method, and the best-fitted b values for detection catalogue and the network catalogue are 1.04±0.01 and 0.79±0.02, respectively.

Fig. 3 The example of template matching method for detecting missing earthquakes (taking the ML2.6 earthquake as an example). Waveform of cross-correlation coefficient (a) after shifting and averaging, blue dots line is the threshold 0.75, red dot is one of the two missing events. (b) shows the waveforms of detected missing earthquake event, red color denotes the template P-wave and blue color corresponding to the S-wave, red dashed line shows the earthquake origin time, the detected event is ML2.6, 07:35:57, 2017-04-16

Fig. 4 Comparison of the frequency and b-value of the detection catalogue and the network catalogue The triangles indicate the detection catalogue, the gray circles represent the network catalogue; Two dashed lines mark the vertex of non-cumulative frequency curves, which represents the magnitude of completeness(Mc). The Mc for detection catalogue and the network catalogue are 0.5 and 1.0, respectively. The b values fitted by the maximum likelihood method are 1.04 and 0.79, respectively
3.2 Characteristics of the Earthquake Swarm

Based on the detection catalogue, we analyze the temporal evolution of the seismicity of Chang Island earthquake swarm. Judging from the variation of cumulative frequency and magnitude with time (Fig. 5), we divide the Chang Island earthquake swarm into four periods: in the first stage, the seismic activity with low magnitudes started to increase from February 13, 2017, the activity continued until the occurrence of the largest ML4.5 earthquake (March 3, 2017). During this time, the network catalog revealed a relative quiet scenario (black point in Fig. 5); the second stage is mainly concentrated between March and April, 2017. After the ML4.5 earthquake on March 3, the seismic activity increased obviously with higher intensity and frequency. However, the seismic activity gradually weakened after reaching the highest value in April, 2017; the third stage occurred in May, 2017, and a new round of enhancement occurred in the earthquake activity from May 2 to 4, and then gradually weakened. The fourth stage is from June to August, 2017. Similar to the third phase, the seismic activity appeared again in the process of enhancement-attenuation, which slowly attenuated after a short time increase in early June. The detection catalogue and the network catalogue have the same seismic activity period, although the number of earthquakes in the low-magnitude section is different. According to the development process of the earthquake swarm, each stage can be further divided into several sub-stages (Fig. 5). In the following, the ETAS model will be used to model these different stages separately; we will obtain the variations of each model parameter and discuss the influence of external force (fluid) action and self-excitation on the development of the earthquake swarm.

Fig. 5 The magnitude-time relation of the network catalogue and the detection catalogue and the cumulative frequency curve of the detection catalogue The left ordinate is the magnitude, the right ordinate is the cumulative frequency, the gray dot indicates the detection catalogue, and the black dot indicates the network catalogue. The long red dotted line divides the earthquake sequence into four stages according to the cumulative frequency. The short red dotted line divides each stage into several sub-stages

The cut-off magnitude has great influence in ETAS modelling (Jiang Haikun et al., 2012; Jiang Changsheng et al., 2013c) and is usually set greater than the magnitude of completeness. To model the seismic activities of all stages, we need to estimate the Mc of the corresponding stage. We use the maximum curvaturemethod and 'goodness of fit' test method to estimate the change of complete magnitude (Mc) with time (Fig. 6). From the results of the optimal solution (the black solid line in Fig. 6), it can be seen that the magnitude of completeness (Mc) is basically constant in the process of swarm development, all are no greater than 0.5. Therefore, in the subsequent parameter inversion of ETAS model, we can set the cut-off magnitude as 0.5.

Fig. 6 The magnitude-sequence map of magnitude of completeness (Mc) The horizontal coordinate represents the earthquake sequence numbe; the gray circle represents earthquake of the detection catalogue; the red line, blue line and cyan dot line represents the magnitude of completeness (Mc) solved by the maximum curvature method, 90% and 95% reliability solved by goodness-of-fit method, respectively. The black solid line represents the optimal solution of the three methods. The window lengths of these methods are all 500 earthquake events
3.3 Inversion of ETAS Model Parameters

According to the stages defined in Fig. 5, the earthquake catalogues of each stage ML≥0.5 are taken, and the respective ETAS model parameters are estimated (Table 1). We select the first and second stages as examples for detailed analysis (Fig. 7).

Table 1 ETAS model parameters at different stages of Chang Island swarm

Fig. 7 Results of ETAS model at different stages. Fig. 7(e) corresponds to the ETAS model parameters of the first stage (Ⅰ) in Fig. 5 and the five sub-phases (Ⅱ1-Ⅱ5) of the second stage, where the horizontal ordinate is the transformed 17:00:48

For the first stage, the ratio of random components is 31.9%, indicating that 31.9% of the seismic activity is triggered by external forces, and 68.1% of the earthquakes are generated by self-excitation of aftershocks. The a-value is 0.608, which is consistent with the activity characteristics of the earthquake cluster (0.35-0.85) (Ogata Y., 1992) and lower than the general tectonic earthquake (2.0). A small a-value indicates that the aftershock excitation ability is not strong. Meanwhile, low α-value means low dependence on earthquake magnitude, therefore, the stress variation caused by the earthquake has only a slight change (Lei Xinglin et al., 2008a). The p-value is larger than the typical tectonic earthquakes (around 1.0), indicating that the sequence decays rapidly. The fitting between theoretical and observed values of the model is poor, and the observed earthquake occurrence rate (blue solid line in Fig. 7(a)) is significantly higher than that of the theoretical curve (grey point line in Fig. 7(a)). This indicates that the self-triggering of Omori's law has only a slight promoting effect, and the proportion of seismic activity caused by external force may be underestimated (Lei Xinglin et al., 2008a).

The second stage corresponds to the period when the seismic activity rate gradually increases and then decreases. According to the activity rate, we divided stage Ⅱ into five sub-phases (Ⅱ1-Ⅱ5, Fig. 5), of which Ⅱ2 has the highest seismic activity rate. The observations and theoretical values of these five sub-phases are well fitted, and the proportion of external force (fluid) action becomes significantly larger and gradually increases, from 46.9% in Ⅱ1 to 63.5% in Ⅱ3, indicating the presence of active external force (fluid). The number of earthquakes triggered by fluids increases, and even more than half of the earthquakes are caused by external forces, instead of being influenced by "aftershocks triggering aftershocks" scenario in Omori's law. At the same time, the corresponding α values are higher than that in the first stage, suggesting that the excitation ability of aftershocks is enhanced. The p-value remains high from stage Ⅰ to Ⅱ3, but decreases significantly in stage Ⅱ4and Ⅱ5, indicating that the decay rate of the sequence decreases.

We calculate the b-value of each stage using maximum likelihood method, and analyze the evolution of the b-value of Chang Island earthquake swarm as well as Rb, α, p and μ from ETAS model in the whole process of swarm development (Fig. 8). The results show that, except for the individual activity stages, the Rb, μ and α are stable and positively correlated. All the values increase at first and then decrease (Fig. 8(a)-(c)), the trend of b-value and p-value are not obviously synchronized, but their value are relative high at early stages and low in later stages. (Fig. 8(d)-(e)). The μ-value is the background seismicity rate. Larger μ value represents stronger external force, and Rb is defined by the ratio of "triggered" earthquakes and "natural" earthquakes unrelated to the external force, which also indicates the triggering intensity (Formula 4). It can be seen from Rb values of different stages that the fluid plays a role in the development of the earthquake swarm to some extent, the effect of fluid increase at first and then weakened, indicating that the extent of external triggering intensity varies at different times.

Fig. 8 Temporal evolution diagrams of ETAS model parameters and b values, in which horizontal ordinate represents 11 stages (including sub-stages) as shown in Fig. 5 and Table 1. Longitudinal coordinates represent changes of parameters
4 DISCUSSION

We use the ratio of earthquakes induced by external force (fluid) and self-excited earthquakes Rb (Formula (4)) to investigate the influence of external force on the swarm activity. Firstly (stage Ⅰ), the proportion of earthquake triggered by external force is not high (31.9%), then the triggering ratio increases gradually until stage Ⅱ3, reaching the maximum (63.5%). Subsequently, the triggering effect of fluids i weakened obvioously.

The b-value usually reflects the stress level of the study area: a low b-value reflects high stress state, while a high b-value corresponds to a low stress state. The high b-value in the early stage of the earthquake swarm indicates that the stress at the beginning of the earthquake swarm is not high (Fig. 8(d)), but the largest ML4.5 earthquake in the swarm occurred, which could be attributed to fluid triggering. Rb indicates the external force (fluid) effect and is generally thought to be related to the stress state of the source region (Peng Yajun et al., 2012; Jia Ke et al., 2014, 2018). In the initial stage of the earthquake swarm, the b-value is high, while the Rb value is relatively low, and the actual earthquake occurrence rate at this stage may be significantly higher than the fitted theoretical curve (Fig. 7(a)) and the external force may be underestimated (Lei Xinglin et al., 2008a). The b-value in other stages is also consistent with the process of increasing at first and then decreasing. Besides, the change of external force is also consistent. When the proportion of fluid action effect is large, the b-value stays in a high level, indicating the fluid reduces the effective normal stress of the fault. Fluids greatly reduce the strength imposed by the fault hence even relatively low local stress can cause earthquakes.

The process of increase and decrease of the earthquake proportion of external force (fluid) triggering may be related to the fluid infiltration process. With the penetration of the fluid, the pore pressure gradually increases, the strength of the fault (fracture) decreases, leading to a larger possibility of earthquake occurrence. Fluid effect at the beginning of the second stage is obviously stronger than that in the firststage; however, after a certain period of time, the fluids tends to be saturated in the pore, and the change of pore pressure decreases gradually, which may weaken the triggering effect of external factors of fluid, and the ratio coefficient Rb decreases from 63.5% to less than 30%. Seismic activities in the third and fourth stages are mainly caused by the self-excitation of earthquakes.

In each stage of swarm development, Rb is positively correlated with μ value, indicating the strength of external force. Larger value corresponds to stronger fluid effect. This is also consistent with the trend of earthquake self-excitation (α). The parameters vary in different stages, however, the fluid interaction may generally have a greater impact on the early stage of the swarm, and the self-excitation of the earthquake plays a more important role in the whole process of swarm development. In the initial stage of the swarm, the p value has been at a high level, result that indicates the swarm attenuates rapidly, suggesting that it may be related to the larger proportion of microearthquakes in the swarm. Alternatively, the p value is positively correlated with the focal depth and the temperature at the focal depth (Creamer F. H. et al., 1993; Nyffenegger P. et al., 2000; Song Jin et al., 2009). Seismic sequences with larger source depths attenuate faster, stress relaxation in regions with higher crustal temperature is faster, and aftershock activity attenuates relatively faster. The focal depth of Chang Island swarm is about 11km, which is larger than that of earthquakes induced by reservoirs and water injection(Jiang Haikun et al., 2012; Lei Xinglin et al., 2008b); Meanwhile, whether the source of fluids is related to high-pressure crustal melting materials in the study area (Liu Lihua et al., 2015) remains questionable. Both of these factors may lead to the accelerated attenuation of earthquake swarms. Detailed changes of background earthquake occurrence rate and p-value can be observed, and subsequent discussions can be completed in the way of iteration and continuous sliding window proposed by Peng Yajun et al. (2012).

5 CONCLUSIONS

In this paper, the GPU-acceleration-based template matching method is used to detect micro-seismic event from the continuous seismic waveforms of the Chang Island earthquake swarm from February 9 to August 20, 2017. A total of 15, 286 earthquake events are identified, of which 302 are self-detection of template events and 14, 984 are missing earthquake events. Events in detected catalogue is more than 6 times of that of the network catalogue, and the magnitude of completeness of the earthquake swarm reduces from 1.0 to 0.5.

Based on the detected earthquake catalogue, the ETAS model is used to analyze the characteristics of the Chang Island earthquake swarm. It is preliminarily inferred that fluid triggering plays a role in the development of the Chang Island swarm to some extent, and has a greater impact in the early stage of the swarm, while the self-excitation of the swarm also plays an important role in the whole development of the swarm. In this paper, the seismogenic mechanism of Chang Island swarm is explored only from the statistical point of view. As for the detailed mechanism and further analyzation of fluid triggering and the change of driving forces in the physical process of swarm development, more methods are needed.

ACKNOWLEDGEMENT

Thanks to the two reviewers for their valuable comments. In this paper, the template matching algorithm is developed by Professor Peng Zhigang, and the ETAS model is calculated by the GeoTaos software developed by Professor Lei Xinglin.

REFERENCES
Akaike H.Akaike H. A new look at the statistical model identification[J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723. DOI:10.1109/TAC.1974.1100705
Creamer F.H., Kisslinger C.Creamer F.H., Kisslinger C. The relation between temperature and the Omori decay parameter for aftershock sequences nearJapan[J]. EOS, 1993, 74(S43): 417.
Gibbons S.J., Ringdal F.Gibbons S.J., Ringdal F. The detection of low magnitude seismic events using array-based waveform correlation[J]. Geophysical Journal International, 2006, 165(1): 149-166. DOI:10.1111/j.1365-246X.2006.02865.x
Hainzl S., Ogata Y.Hainzl S., Ogata Y. Detecting fluid signals in seismicity data through statistical earthquake modeling[J]. Journal of Geophysical Research:Solid Earth, 2005, 110(B5): B05S07. DOI:10.1029/2004JB003247
Helmstetter A., Sornette D., Grasso J.R.Helmstetter A., Sornette D., Grasso J.R. Mainshocks are aftershocks of conditional foreshocks:How do foreshock statistical properties emerge from aftershock laws[J]. Journal of Geophysical Research:Solid Earth, 2003, 108(B1): 2046. DOI:10.1029/2002JB001991
Hou Jinxin, Wang BaoshanHou Jinxin, Wang Baoshan. Temporal evolution of seismicity before and after the 2014 Ludian Ms6.5 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(4): 1446-1456 (in Chinese with English abstract).
Huang Yilei, Zhou Shiyong, Zhuang JiancangHuang Yilei, Zhou Shiyong, Zhuang Jiancang. Numerical tests on catalog-based methods to estimate magnitude of completeness[J]. Chinese Journal of Geophysics, 2016, 59(4): 1350-1358 (in Chinese with English abstract).
Jia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng, Guo Yicun, Gao Zhaohui, Gao SheshengJia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng, Guo Yicun, Gao Zhaohui, Gao Shesheng. Did the 2008 MW7.9 Wenchuan earthquake trigger the occurrence of the 2017 MW6.5 Jiuzhaigou earthquake in Sichuan, China?[J]. Journal of Geophysical Research:Solid Earth, 2018, 123(4): 2965-2983. DOI:10.1002/2017JB015165
Jia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang ChangshengJia Ke, Zhou Shiyong, Zhuang Jiancang, Jiang Changsheng. Possibility of the Independence between the 2013 Lushan Earthquake and the 2008 Wenchuan Earthquake on Longmen Shan Fault, Sichuan, China[J]. Seismological Research Letters, 2014, 85(1): 60-67. DOI:10.1785/0220130115
Jiang Changsheng, Zhuang Jiancang, Long Feng, Han Libo, Guo LujieJiang Changsheng, Zhuang Jiancang, Long Feng, Han Libo, Guo Lujie. Statistical analysis of ETAS parameters in the early stage of the 2013 Lushan MS7.0 earthquake sequence[J]. Acta Seismologica Sinica, 2013a, 35(5): 661-669 (in Chinese with English abstract).
Jiang Changsheng, Wu Zhongliang, Zhuang JiancangJiang Changsheng, Wu Zhongliang, Zhuang Jiancang. ETAS model applied to the Earthquake-Sequence Association (ESA) problem:the Tangshan sequence[J]. Chinese Journal of Geophysics, 2013b, 56(9): 2971-2981 (in Chinese with English abstract).
Jiang Changsheng, Wu Zhongliang, Han Libo, Guo LujieJiang Changsheng, Wu Zhongliang, Han Libo, Guo Lujie. Effect of cutoff magnitude Mc of earthquake catalogues on the early estimation of earthquake sequence parameters with implication for the probabilistic forecast of aftershocks:the 2013 Minxian-Zhangxian, Gansu, MS6.6 earthquake sequence[J]. Chinese Journal of Geophysics, 2013c, 56(12): 4048-4057 (in Chinese with English abstract).
Jiang Changsheng, Wu Zhongliang, Yin Fengling, Guo Lujie, Bi Jinmeng, Wang YawenJiang Changsheng, Wu Zhongliang, Yin Fengling, Guo Lujie, Bi Jinmeng, Wang Yawen. Stability of early-estimation sequence parameters for continuous forecast of the aftershock rate:a case study of the 2014 Ludian, Yunnan MS6.5 earthquake[J]. Chinese Journal of Geophysics, 2015, 58(11): 4163-4173 (in Chinese with English abstract).
Jiang Changsheng, Zhuang Jiancang, Wu Zhongliang, Bi JinmengJiang Changsheng, Zhuang Jiancang, Wu Zhongliang, Bi Jinmeng. Application and comparison of two short-term probabilistic forecasting models for the 2017 Jiuzhaigou, Sichuan, MS7.0 earthquake[J]. Chinese Journal of Geophysics, 2017, 60(10): 4132-4144 (in Chinese with English abstract).
Jiang Haikun, Zheng Jianchang, Wu Qiong, Qu Yanjun, Li YongliJiang Haikun, Zheng Jianchang, Wu Qiong, Qu Yanjun, Li Yongli. Earlier statistical features of ETAS model parameters and their seismological meanings[J]. Chinese Journal of Geophysics, 2007, 50(6): 1778-1786 (in Chinese with English abstract).
Jiang Haikun, Yang Maling, Sun Xuejun, Lü Jian, Yan Chunheng, Wu Qiong, Song Jin, Zhao Yong, Huang Guohua, Zhang Hua, Yao Hong, Mu Jianying, Li Jin, Qu JunhaoJiang Haikun, Yang Maling, Sun Xuejun, Lü Jian, Yan Chunheng, Wu Qiong, Song Jin, Zhao Yong, Huang Guohua, Zhang Hua, Yao Hong, Mu Jianying, Li Jin, Qu Junhao. A typical example of locally triggered seismicity in the boundary area of Lingyun and Fengshan following the large rainfall event of June 2010[J]. Chinese Journal of Geophysics, 2011, 54(10): 2606-2619 (in Chinese with English abstract).
Jiang Haikun, Song Jin, Wu Qiong, Li Jin, Qu JunhaoJiang Haikun, Song Jin, Wu Qiong, Li Jin, Qu Junhao. Quantitative investigation of fluid triggering on seismicity in the Three Gorge Reservoir area based on ETAS model[J]. Chinese Journal of Geophysics, 2012, 55(7): 2341-2352 (in Chinese with English abstract).
Lei Xinglin, Ma Shengli, Wen Xueze, Su Jinrong, Du FangLei Xinglin, Ma Shengli, Wen Xueze, Su Jinrong, Du Fang. Integrated analysis of stress and regional seismicity by surface loading-a case study of Zipingpu reservoir[J]. Seismology and Geology, 2008a, 30(4): 1046-1064 (in Chinese with English abstract).
Lei Xinglin, Yu Guozheng, Ma Shengli, Wen Xueze, Wang QiangLei Xinglin, Yu Guozheng, Ma Shengli, Wen Xueze, Wang Qiang. Earthquakes induced by water injection at ~3km depth within the Rongchang gas field, Chongqing, China[J]. Journal of Geophysical Research:Solid Earth, 2008b, 113(B10): B10310. DOI:10.1029/2008JB005604
Lei Xinglin, Huang Dongjian, Su Jinrong, Jiang Guomao, Wang Xiaolong, Wang Hui, Guo Xin, Fu HongLei Xinglin, Huang Dongjian, Su Jinrong, Jiang Guomao, Wang Xiaolong, Wang Hui, Guo Xin, Fu Hong. Fault reactivation and earthquakes with magnitudes of up to Mw4.7 induced by shale-gas hydraulic fracturing in Sichuan Basin, China[J]. Scientific Reports, 2017, 7(1): 7971. DOI:10.1038/s41598-017-08557-y
Li Zhichao, Huang QinghuaLi Zhichao, Huang Qinghua. Assessment of detectability of the Capital-circle Seismic Network by using the probability-based magnitude of completeness (PMC) method[J]. Chinese Journal of Geophysics, 2014, 57(8): 2584-2593 (in Chinese with English abstract).
Liu Lihua, Hao Tianyao, Lü Chuanchuan, You Qingyu, Pan Jun, Wang Fuyun, Xu Ya, Zhao Chunlei, Zhang JianshiLiu Lihua, Hao Tianyao, Lü Chuanchuan, You Qingyu, Pan Jun, Wang Fuyun, Xu Ya, Zhao Chunlei, Zhang Jianshi. Crustal structure of Bohai Sea and adjacent area (North China) from two onshore-offshore wide-angle seismic survey lines[J]. Journal of Asian Earth Sciences, 2015, 98: 457-469. DOI:10.1016/j.jseaes.2014.11.034
Long Feng, Du Fang, Ruan Xiang, Deng Yuanqing, Zhang TiebaoLong Feng, Du Fang, Ruan Xiang, Deng Yuanqing, Zhang Tiebao. Water injection triggered earthquakes in the Zigong mineral wells in ETAS model[J]. Earthquake Research in China, 2010, 26(2): 164-171 (in Chinese with English abstract).
Ma Tengfei. "Repeating Earthquakes" Extracted from Continuous Waveforms: Analysis of the Wenchuan Sequence in Connection to the WIFSD Drilling Site[D]. Beijing: Institute of Geophysics, China Earthquake Administration, 2015 (in Chinese).
Meng Xiaofeng, Yu Xiao, Peng Zhigang, Hong BoMeng Xiaofeng, Yu Xiao, Peng Zhigang, Hong Bo. Detecting earthquakes around Salton Sea following the 2010 MW7.2 El Mayor-Cucapah earthquake using GPU parallel computing[J]. Procedia Computer Science, 2012, 9: 937-946. DOI:10.1016/j.procs.2012.04.100
Narteau C., Byrdina S., Shebalin P., Schorlemmer D.Narteau C., Byrdina S., Shebalin P., Schorlemmer D. Common dependence on stress for the two fundamental laws of statistical seismology[J]. Nature, 2009, 462(7273): 642-646. DOI:10.1038/nature08553
Nyffenegger P., Frohlich C.Nyffenegger P., Frohlich C. Aftershock occurrence rate decay properties for intermediate and deep earthquake sequences[J]. Geophysical Research Letters, 2000, 27(8): 1215-1218. DOI:10.1029/1999GL010371
Ogata Y.Ogata Y. Statistical models for earthquake occurrences and residual analysis for point processes[J]. Journal of the American Statistical Association, 1988, 83(401): 9-27. DOI:10.1080/01621459.1988.10478560
Ogata Y.Ogata Y. Detection of precursory relative quiescence before great earthquakes through a statistical model[J]. Journal of Geophysical Research:Solid Earth, 1992, 97(B13): 19845-19871. DOI:10.1029/92JB00708
Ogata Y.Ogata Y. Increased probability of large earthquakes near aftershock regions with relative quiescence[J]. Journal of Geophysical Research:Solid Earth, 2001, 106(B5): 8729-8744. DOI:10.1029/2000JB900400
Omori F.Omori F. On the aftershocks of earthquake[J]. Journal of the College of Science, Imperial University of Tokyo, 1894, 7: 11-200.
Peng Yajun, Zhou Shiyong, Zhuang Jiancang, Shi JiaPeng Yajun, Zhou Shiyong, Zhuang Jiancang, Shi Jia. An approach to detect the abnormal seismicity increase in Southwestern China triggered co-seismically by 2004 Sumatra MW9.2 earthquake[J]. Geophysical Journal International, 2012, 189(3): 1734-1740. DOI:10.1111/j.1365-246X.2012.05456.x
Peng Zhigang, Zhao PengPeng Zhigang, Zhao Peng. Migration of early aftershocks following the 2004 Parkfield earthquake[J]. Nature Geoscience, 2009, 2(12): 877-881. DOI:10.1038/ngeo697
Schaff D.P.Schaff D.P. Semiempirical statistics of correlation-detector performance[J]. Bulletin of the Seismological Society of America, 2008, 98(3): 1495-1507. DOI:10.1785/0120060263
Shelly D.R., Beroza G.C., Ide S.Shelly D.R., Beroza G.C., Ide S. Non-volcanic tremor and low-frequency earthquake swarms[J]. Nature, 2007, 446(7133): 305-307. DOI:10.1038/nature05666
Song Jin, Jiang HaikunSong Jin, Jiang Haikun. A review on decay and generation of aftershock activity[J]. Seismology and Geology, 2009, 31(3): 559-571 (in Chinese with English abstract).
Su LuanshengSu Luansheng. Relationship between the fractal structure variation in Changdao area and the arealsei smicity and short-term and impending precursor feature[J]. Journal of Seismological Research, 1993, 16(1): 33-40 (in Chinese with English abstract).
Tan Yipei, Deng Li, Cao Jingquan, Shan LianjunTan Yipei, Deng Li, Cao Jingquan, Shan Lianjun. Seismological mechanism analysis of 2015 Luanxian swarm, Hebei province[J]. Chinese Journal of Geophysics, 2016, 59(11): 4113-4125 (in Chinese with English abstract).
Utsu T.Utsu T. A statistical study on the occurrence of aftershocks[J]. Geophysical Magazine, 1961, 30: 526-605.
Wang Peng, Hou Jinxin, Wu PengWang Peng, Hou Jinxin, Wu Peng. Temporal evolution of the seismicity of the 2017 Jiuzhaigou MS7.0 earthquake sequence[J]. Earthquake Research in China, 2017, 33(4): 453-462 (in Chinese with English abstract).
Wang Zhicai, Deng Qidong, Chao Hongtai, Du Xiansong, Shi Ronghui, Sun Zhaomin, Xiao Lanxi, Min Wei, Ling HongWang Zhicai, Deng Qidong, Chao Hongtai, Du Xiansong, Shi Ronghui, Sun Zhaomin, Xiao Lanxi, Min Wei, Ling Hong. Shallow-depth sonic reflection profiling studies on the active Penglai-Weihai fault zone offshore of the northern Shandong peninsula[J]. Chinese Journal of Geophysics, 2006, 49(4): 1092-1101 (in Chinese with English abstract).
Wiemer S., Wyss M.Wiemer S., Wyss M. Minimum magnitude of completeness in earthquake catalogs:Examples from Alaska, the Western United States, and Japan[J]. Bulletin of the Seismological Society of America, 2000, 90(4): 859-869. DOI:10.1785/0119990114
Xu Jie, Song Changqing, Chu QuanzhiXu Jie, Song Changqing, Chu Quanzhi. Preliminary study on the Seismotectonic characters of the Zhangjiakou Penglai fault zone[J]. Seismology and Geology, 1998, 20(2): 146-154 (in Chinese with English abstract).
Yang Hongfeng, Zhu Lupei, Chu RishengYang Hongfeng, Zhu Lupei, Chu Risheng. Fault-plane determination of the 18 April 2008 Mount Carmel, Illinois, earthquake by detecting and relocating aftershocks[J]. Bulletin of the Seismological Society of America, 2009, 99(6): 3413-3420. DOI:10.1785/0120090038
Zhang Miao, Wen LianxingZhang Miao, Wen Lianxing. Seismological evidence for a low-yield nuclear test on 12 May 2010 in North Korea[J]. Seismological Research Letters, 2015, 86(1): 138-145. DOI:10.1785/02201401170
Zheng Jianchang, Xu Changping, Li Dongmei, Zhou CuiyingZheng Jianchang, Xu Changping, Li Dongmei, Zhou Cuiying. Recent change of stress field in North China derived from station composite focal mechanisms[J]. Recent Developments in World Seismology, 2018(4): 48-57 (in Chinese with English abstract).
Zhuang JiancangZhuang Jiancang. Statistical modeling of seismicity patterns before and after the 1990 Oct 5 Cape Palliser earthquake, New Zealand[J]. New Zealand Journal of Geology and Geophysics, 2000, 43(3): 447-460. DOI:10.1080/00288306.2000.9514901
Zhuang Jiancang, Ogata Y., Wang TingZhuang Jiancang, Ogata Y., Wang Ting. Data completeness of the Kumamoto earthquake sequence in the JMA catalog and its influence on the estimation of the ETAS parameters[J]. Earth, Planets and Space, 2017, 69(1): 36.